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Abstract 

Mathematical  modeling  and  atomic  force  microscopy  were  used  to  optimize  a 
multi-layer  reed  sensor  for  use  in  a  new  infrared  detection  system.  The  reed  sensor  is 
designed  to  deflect  in  response  to  temperature  variations.  In  this  system,  a  sample 
absorbs  infrared  light  (modulated  at  some  frequency  to  reduce  noise)  and  releases  heat 
into  the  reed.  The  layers  expand  differently,  causing  the  reed  to  deflect.  The  deflection 
is  influenced  by  geometric  parameters,  such  as  length  and  thickness  of  each  layer,  and 
also  by  the  frequency  of  modulation.  Deflection  is  easily  measured  using  an  atomic  force 
microscope.  This  research  applied  a  mathematical  model  to  the  reed  sensor,  explored  the 
effect  of  varying  geometric  parameters,  and  used  experimental  testing  to  evaluate  the 
model.  Several  reeds  were  tested  to  validate  the  model,  which  predicted  sensitivity, 
maximum  deflection,  resonance  frequency,  and  signal  to  noise  ratio.  Varying 
concentrations  of  absorbing  chemical  were  deposited  in  the  reed,  to  cause  a  change  in 
deflection. 

A  cantilever  reed  having  three  layers  was  modeled  mathematically  using  a 
recently  developed  thermal  model  and  a  mechanical  model  based  on  two-layer  theory. 

The  predictions  mirror  theory-based  expectations  in  that  for  each  reed  there  was  an 
optimum  modulating  frequency.  Also,  the  model  was  able  to  predict  optimum  thickness 
ratios,  and  was  used  to  gain  insight  into  the  phase  shift  behavior  of  reed  deflection. 

This  model  behavior  is  in  line  with  laboratory  experiments,  with  some  exceptions. 
The  model  correctly  predicts  temperature  behavior,  but  currently  does  not  provide 
accurate  quantitative  results.  The  model  is  found  lacking  for  short,  squarer  reeds  not 
having  the  assumed  aspect  ratio  of  10.  This  type  of  reed  oscillated  in  higher  order  modes 
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than  were  considered  in  the  mathematical  model.  To  accurately  represent  these,  the 
thermal  and  mechanical  behavior  must  be  derived  using  plate  theory  rather  than  beam 
theory.  The  model  is  also  inadequate  in  describing  the  effects  of  thickening  the  top  layer, 
probably  because  it  assumes  uniform  energy  deposition  throughout  the  layer.  Finally, 
the  deflection  signal  was  observed  to  decrease  with  increased  concentration  of  sample 
deposited. 


Sponsor:  US  Army  CBDCOM 
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OPTIMIZATION  OF  A  MULTILAYER  PHOTOTHERMAL 
SENSOR  FOR  INFRARED  SPECTROSCOPY 

1.  Introduction 


Background 

Several  branches  of  the  military  have  joined  in  an  effort  to  improve  their  water 
contamination  detection  systems.  Some  important  system  requirements  are  summarized 
below  [10:  1-3]: 

-Detect  and  quantify  specified  chemical  and  biological  agents 

-Achieve  Parts-Per-Billion  (PPB)  Detection  level 

-Field  operability  by  military  personnel 

-Be  portable  by  one  individual. 

This  system  must  be  rugged,  small,  and  have  rapid  measurement  capability. 
Infrared  detectors  have  good  selectivity  for  measuring  a  wide  variety  of  agents.  Because 
they  are  based  on  measuring  very  small  differences  in  large  signals,  current  infrared 
detectors  are  too  noisy  at  room  temperature  to  be  useful  in  this  military  application.  AFIT 
proposes  to  meet  this  challenge  via  a  new  photothermal  infrared  detection  system. 

The  principle  behind  this  new  system  is  simple.  A  sample  is  captured  in 
nanopores  in  magnesium  aluminate,  which  is  the  top  layer  of  a  cantilever  reed  (refer  to 
Figure  1  below).  To  search  for  a  particular  contaminant,  light  of  an  appropriate 
wavelength  is  directed  onto  the  reed.  Nanoporous  magnesium  aluminate,  also  known  as 
sol-gel,  was  chosen  as  the  top  layer  because  it  is  transparent  in  the  IR  and  near  IR.  The 
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chemicals  of  interest  to  the  military  are  typically  responsive  to  wavelengths  in  IR  and 
near  IR  light.  If  the  contaminant  is  present,  light  will  be  absorbed  and  re-released  as  heat. 
Heat  flows  to  the  other  layers  of  the  cantilever  reed,  which  expand  differently,  according 
to  the  difference  in  their  thermal  expansion  coefficients.  Because  the  layers  are  tightly 
bound  together  at  the  interfaces,  the  net  effect  is  the  bending  of  the  reed.  Deflection, 
then,  should  be  a  direct  measure  of  how  much  contaminant  is  present  in  the  sol-gel  layer. 

The  novel  part  of  the  system  is  the  reed  upon  which  those  contaminants  are 
deposited.  The  reed  is  made  of  several  layers,  beginning  with  a  nanoporous  magnesium 
aluminate  (sol-gel,  MgAl2C>4)  layer,  typically  one  micron  in  thickness.  Next  is  a  thin 
(about  200  Angstroms)  layer  of  gold,  followed  by  layers  of  silicon,  typically  eight 
microns,  and  aluminum,  typically  half  a  micron.  The  sol-gel  layer  entraps  molecules  of 
the  sample  to  be  examined;  the  gold  layer  conducts  heat  to  the  remainder  of  the  reed,  and 
effectively  isolates  the  silicon  and  aluminum  layers  from  the  stiff  sol-gel  layer. 

Deflection  is  dominated  by  the  expansion  of  the  bottom  two  layers,  with  aluminum 
expanding  about  ten  times  more  than  silicon.  Deflection  measurement  response  to 
wavelength  of  light  provides  a  quantitative  measurement  of  the  quantity  of  pollutant 
trapped.  The  figure  below  (Figure  1.1)  shows  this  multi-layer  reed  as  used  in  the 
laboratory.  The  sample  is  located  in  the  sol-gel  layer.  The  reed  is  held  to  its  silicon  base 
by  a  clamp  and  is  covered  by  a  reflector,  which  forces  the  incoming  light  to  make 
repeated  passes  over  the  reed  sensor. 
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Figure  1.1:  Cantilever  Reed  Design  (not  to  scale) 

The  sol-gel  material  was  chosen  for  the  first  layer  because  it  was  available  in 
nanoporous  form  (which  easily  traps  the  sample  molecules)  and  because  it  is  relatively 
transparent  to  infrared  light  [11:  38].  Silicon  also  absorbs  infrared  light,  as  will  be 
discussed  later. 

There  are  several  systems  available  for  measurement  of  the  cantilever  reed 
deflection.  In-house  testing  using  a  fiber  interferometer  was  successful,  but  complicated 
compared  to  measurements  with  an  atomic  force  microscope  (AFM).  In  this  project, 
deflection  was  measured  using  the  AFM,  which  employs  a  photodiode  and  laser  beam 
(combined  to  produce  a  position-sensitive  detector),  focused  on  the  100-micron  AFM 
cantilever  tip  in  contact  with  the  end  of  the  reed  sensor.  This  system  is  small,  does  not 
need  the  cryogenics  required  by  current  infrared  detectors,  and  is  expected  to  detect 
samples  at  parts-per-billion  concentrations.  Thus,  the  system  as  it  currently  stands  meets 
all  of  the  system  requirements  listed  earlier,  with  the  exception  of  the  parts-per-billion 
level  of  detection.  This  requirement  should  be  met  once  the  system  has  been  optimized. 
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Essential  in  producing  the  best  possible  detector  is  the  optimum  reed  design.  The 
best  reed  is  as  thin  as  possible,  while  maintaining  rigidity.  Also,  it  must  be  long  and 
wide  enough  to  capture  a  sufficient  amount  of  sample,  but  not  so  long  as  to  be  overly 
sensitive  to  low  frequency  noise.  Geometric  parameters  such  as  thickness,  length  and 
width  of  each  layer  will  be  varied  and  the  performance  of  each  reed  will  be  evaluated 
using  models  and  laboratory  testing. 

Problem  Statement 

This  research  is  concerned  with  optimizing  a  cantilever  reed  for  use  in  this 
photothermal  infrared  measurement  system.  Mathematical  modeling  and  laboratory 
testing  of  these  reeds  is  used  to  gain  an  understanding  of  how  a  change  in  the  cantilever 
reed  geometry  affects  sensitivity.  The  goal  is  to  increase  sensitivity. 

Research  Focus 

The  cantilever  reed  was  modeled  mathematically  and  the  deflection  response  to 
changes  in  such  parameters  as  length  and  thickness  of  layers  were  calculated.  These 
models,  which  include  a  thermal  model  of  the  three  important  layers  and  a  mechanical 
model  based  on  two-layer  theory,  predicted  sensitivity,  maximum  deflection,  resonance 
frequency,  and  signal  to  noise  ratio,  all  of  which  must  be  balanced  when  determining  the 
optimum  reed  parameters.  Several  reeds,  of  varying  dimensions,  were  constructed  and 
tested  to  validate  the  mathematical  model. 

Research  Concerns/Hypotheses 
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The  first  step  in  increasing  sensitivity  of  the  reed  is  minimizing  noise  sources. 
Absorption  of  the  IR  light  by  the  silicon  layer,  rather  than  by  the  sample  trapped  within 
the  sol-gel  layer,  produces  a  background  signal  many  times  larger  than  the  deflection 
signals  due  to  the  sample.  Silicon  shows  strong  absorption  in  the  IR,  with  a  prominent 
band  at  1.2  microns. 

A  silicon  filter  placed  between  the  light  source  and  the  three-layer  reed  can  reduce 
this  problem,  by  removing  the  wavelengths  of  light  absorbed  by  silicon.  This  would 
reduce  the  overall  deflection,  because  less  light  is  hitting  the  reed,  but  would  increase  the 
percentage  of  deflection  associated  with  the  sample  itself. 

One  proposal  to  prevent  the  IR  light  from  exciting  the  silicon  was  to  use  a  layer 
of  gold  between  the  sol-gel  and  the  silicon  layers.  This  gold  layer  would  reflect  the  IR 
beam  back  through  the  sol-gel.  This  turned  out  to  require  too  thick  a  layer,  as  thick  as  the 
wavelength  of  light  used.  The  presence  of  the  gold,  however,  helped  eliminate  the 
following  mechanical  problem.  The  sol-gel  material  has  twice  the  thermal  expansion  of 
silicon.  Table  A.l  of  the  Appendix  A  presents  material  properties  for  each  of  the  layers. 
The  sol-gel  is  also  expected  to  heat  up  the  most  because  it  receives  the  incoming  light 
directly  but  loses  heat  slowly  because  of  its  small  thermal  conductivity  and  diffusivity. 

The  sol-gel  will  expand  more  than  the  silicon.  This  would  cause  downward  deflection, 
unless  the  aluminum  absorbs  enough  energy  to  force  the  reed  upward.  This  concern  is 
practically  eliminated  by  the  use  of  the  gold  layer,  which  is  assumed  to  be  ductile  enough 
to  mechanically  isolate  the  stiff  sol-gel  layer  from  the  silicon/aluminum  combination. 

This  assumption  of  mechanical  isolation  of  the  sol-gel,  due  to  the  presence  of  the 
gold  layer,  also  simplifies  modeling.  It  is  thin  enough  and  conductive  enough  that  it  need 
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not  be  included  in  thermal  modeling  of  the  reed.  Because  gold  is  both  ductile  and 
conductive,  it  essentially  separates  the  sol-gel  layer  from  the  other  layers,  while 
efficiently  conducting  heat  into  the  silicon.  Thus,  two-layer  mechanical  models  may  be 
used,  circumventing  development  of  a  more  complicated  multi-layer  model.  Of  course, 
the  solution  still  requires  a  three-layer  thermal  solution  to  determine  correct  temperature 
distribution. 

Sensitivity  of  the  detector  is  increased  when  operating  at  the  mechanical  resonant 
frequency  of  the  reed.  At  this  frequency,  the  reed  deflects  strongly  with  even  small 
changes  in  temperature.  Operating  near  this  frequency  makes  the  reed  more  sensitive  to 
both  temperature  change  and  noise  [12:  2486].  Comparison  of  signal-to-noise  ratio 
(SNR)  on  and  off  mechanical  resonance  frequency  was  expected  to  favor  operation  at  the 
mechanical  resonance. 

Another  way  to  increase  the  sensitivity  of  the  reed  is  to  increase  its  length. 
Deflection  and  sensitivity  (deflection  divided  by  input  energy)  are  both  directly  related  to 
length  [4:  32].  The  difference  in  length  between  expanded  layers  will  be  increased  for  a 
longer  reed,  and  the  longer  reed  will  have  a  larger  surface  area  to  trap  more  pollutant,  if 
other  factors  are  held  constant.  Increasing  the  length  means  that  each  layer  will  expand 
more,  resulting  in  greater  deflection.  Unfortunately,  lengthening  the  reed  can  lower  the 
resonant  frequency,  making  the  reed  is  extremely  sensitive  to  low  frequency  noise. 
Methodology 

These  effects  are  better  understood  by  using  a  mathematical  model  to  describe 
behavior  of  the  reed.  First,  the  diffusion  of  heat  between  each  of  the  layers  must  be 
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determined.  The  mechanical  portion  of  the  problem  then  uses  the  thermal  solution  to 
determine  deflection  and  phase  lag  information. 

Current  two-layer  reed  thermal  models  are  inadequate  because  they  describe  only 
the  static  case  [2:  3793, 15:  83].  A  dynamic  solution  is  required  to  describe  the  deflection 
amplitude  and  phase  lag  resulting  from  the  reed’s  dynamic  temperature  change. 
Mechanical  isolation  of  the  sol-gel  allows  use  of  the  two-layer  mechanical  equations, 
given  the  correct  three-layer  temperature  distribution. 

A  model  has  recently  been  developed  by  Dr  William  Baker,  Associate  Professor 
of  Mathematics  at  the  Air  Force  Institute  of  Technology,  to  describe  thermal  behavior  of 
a  three-layer  reed.  It  is  presented  in  Chapter  2.  Once  the  model  has  been  used  to  predict 
optimum  dimensions,  several  reeds  near  the  optimum  reed  size  can  be  fabricated  and 
tested. 

Testing.  Error  in  deflection  measurement  was  reduced  using  a  chopper  to 
modulate  the  IR  beam.  A  lock-in  amplifier  provided  phase  and  amplitude  information. 
Later  in  the  project,  faster  digital  filter  methods  were  developed  to  improve  accuracy  of 
deflection  amplitude  measurement.  Amplitude  of  the  output  signal  is  proportional  to  the 
amplitude  of  the  deflection,  as  expected.  Phase  information  is  related  to  the  rate  of  heat 
transfer  through  the  sol-gel  layer.  Phase  lag  arises  from  the  thermal  diffusion  rate  through 
layers  of  the  reed.  The  sol-gel  makes  the  largest  contribution  to  phase  lag,  because  of  its 
small  thermal  diffusion.  Phase  lag  increases  with  increased  modulation  frequency  of  the 
IR  light. 

The  modulation  frequency  was  varied  and  the  corresponding  deflection  and  phase 
lag  responses  were  observed.  Varying  concentrations  of  absorbing  chemical  were 
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deposited  in  the  reed,  to  cause  a  change  in  deflection.  Several  reed  lengths  were  tested, 
with  the  shorter  reeds  expected  to  produce  less  noisy  deflection  data. 

Integration  of  Reed  Into  System.  Originally,  adhesive  was  used  to  attach  the  reed 
to  the  base,  but  the  methanol  solvent  used  to  dissolve  the  absorber  (pyridine)  dissolved 
the  glue  and  caused  the  reeds  to  fall  off.  Currently,  a  clamp  is  used  to  attach  the  reed  to 
its  base.  This  method  is  undesirable  because  the  clamps  themselves  tend  to  shear  the 
reeds  at  the  point  of  contact.  Although  this  loss  can  be  minimized  with  careful  closing  of 
the  clamp,  it  is  better  to  find  a  less  forceful  method  of  keeping  the  reed  on  the  base. 

A  mirror/reflector  can  increase  the  deflection  signal  from  the  reed.  This  mirror, 
typically  placed  above  the  reed,  would  reflect  the  IR  beam  back  toward  the  reed.  The 
light  beam  hits  the  cantilever  reed  many  times,  depositing  more  energy,  and  improving 
sensitivity.  Re-usable  mirrors  of  polished  brass  and  aluminum  were  created  and  tested. 

Silicon  is  particularly  absorptive  in  the  near  IR.  Heating,  due  to  silicon 
absorption,  causes  a  deflection  response  resulting  in  higher  background.  A  silicon 
window  may  solve  this  problem  by  absorbing  light  from  the  prominent  silicon  absorption 
bands  before  the  light  hits  the  reed.  This  could  reduce  the  background  signal  when 
exploring  reed  properties  with  a  broadband  source.  Because  no  mid-IR  tunable  laser  was 
available,  a  near  IR  (from  visible  to  about  8  microns)  broad  band  source  was  used  in  this 
research. 

Assumptions/Limitations 

The  mathematical  model  assumes  that  the  temperature  profile  of  the  reed  does  not 
vary  with  width,  and  that  energy  is  deposited  uniformly  in  the  first  layer.  Convection  and 
radiation  effects  are  ignored,  as  are  noise  sources  other  than  thermal  noise.  The  use  of 
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beam  theory  to  describe  the  reed  sensor  requires  the  aspect  ratio  (length  /  width)  will 
always  be  at  least  ten  and  the  length  to  be  24  times  greater  than  thickness.  The  latter 
condition  is  always  met  within  the  parameter  space  of  this  problem.  This  former 
condition  is  necessary  to  assume  that  deflection  occurs  only  as  a  function  of  distance 
down  the  length  of  the  reed. 

Design  constraints  include  the  area  and  wavelength  range  of  the  light  beam 
impinging  the  reed.  Width  of  the  reed  is  bounded  on  the  lower  end,  by  the  diameter  of 
the  light  beam,  currently  half  of  a  millimeter.  A  laser  source  can  be  focused  down  to  0.01 
mm  or  less.  At  these  small  widths,  production  constraints  define  the  boundaries.  The 
smallest  reed  that  can  be  prepared  using  current  techniques  is  half  a  millimeter  in  width 
and  3  millimeters  in  length. 

Implications 

Optimizing  the  cantilever  reed  is  an  important  step  in  the  process  of  optimizing 
this  photothermal  infrared  detector.  The  results  validate  Baker’s  currently  unpublished 
three-layer  thermal  equations,  presented  in  Appendix  B,  which  can  now  be  combined 
with  other  models  so  that  behavior  of  even  smaller  reed  sensors  can  be  determined.  This 
research  helps  AFIT  achieve  its  goals  by  making  the  detector  more  sensitive. 
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2.  Review 


Historical  Perspective 

This  is  a  new  application  for  the  position  sensitive  detection  implemented  in  the 
atomic  force  microscope  (AFM).  Researchers  are  just  now  beginning  to  adapt  the 
cantilever  of  the  AFM  itself  for  use  in  temperature  detection  systems. 

The  AFM  measures  deflection  by  measuring  the  movement  of  a  cantilever  whose 
small  tip  will  be  placed  on  the  end  of  the  multi-layer  reed.  As  the  reed  deflects,  the  tip  is 
moved  up  and  down.  A  laser  beam  directed  at  the  AFM  cantilever  reflects  to  a 
photodiode,  which  tracks  movement  of  the  laser,  cantilever  and  therefore  the  reed.  AFM 
cantilever  properties  have  been  well  established,  and  many  of  the  resulting  models  apply 
to  other  vibrating  levers  and  so  can  be  used  in  describing  the  sample  laden,  multi-layer 
reed. 

In  1997,  Moulin,  Stephenson  and  Welland  ran  a  three  dimensional  model  of  a 
two-layer  triangular  shaped  AFM  cantilever  in  order  to  quantify  its  mechanical  and 
thermal  properties  when  operating  in  either  air  or  water.  They  were  able  to  accurately 
predict  time  constant  and  deflection  for  bi-layer  reeds  [9:  590],  Conventional  models 
were  using  only  one  of  the  reed’s  resonance  frequencies  and  did  not  account  for  two- 
dimensional  vibration  of  the  reed.  Salapka,  et  al,  used  a  multi-mode  model  to 
characterize  the  noise  present  in  such  measurements,  including  both  additional 
resonances,  and  higher  modes.  These  factors  were  necessary  to  accurately  predict  noise 
in  reeds  of  aspect  ratio  below  eight  [12:  2533].  Walters,  et  al,  explored  the  value  of 
operating  in  the  cantilever’s  mechanical  resonance  region  and  found  an  increase  in  signal 
to  noise  (S/N)  ratio  in  spite  of  a  marked  increase  in  noise  density  [14:  3587].  Lai,  et  al, 
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optimized  a  two-layer  cantilever  for  use  in  a  similar  thermal  sensor.  They  coated  a 
silicon  AFM  cantilever  with  a  second  layer  and  measured  its  deflection  response  to 
change  in  temperature.  By  optimizing  the  thickness  of  the  layers,  they  improved 
sensitivity  by  “40%  over  previous  studies”  [6:  117]. 

One  notable  bit  of  research  is  that  of  J.  Barnes,  et  al  reported  in  1994.  They  used 
a  three-layer  cantilever  in  the  AFM  to  measure  the  presence  of  a  chemical  [3:  79].  The 
multi-layer  cantilever  was  made  from  a  commercially  prepared  two-layer  AFM  cantilever 
(silicon  nitride  and  aluminum)  with  140  nm  latex  spheres  deposited  on  top.  They  were 
able  to  detect  samples  with  this  method,  claiming  a  demonstrated  sensitivity  a  factor  of 
100  higher  than  previous  methods.  Though  Bames  (et  al)  measured  deflection  and 
resonance  response  with  a  third  layer,  they  presented  no  mathematical  description  of  their 
reed. 

Finally,  Burggraf  and  Li  extended  the  idea  to  the  use  of  a  transparent  sol-gel  layer 
on  a  separate  reed  [4:  1],  The  use  of  an  additional  reed  increases  throughput  so  that 
different  samples  can  easily  be  measured  in  succession.  The  MgAl204  sol-gel  layer 
appears  to  be  a  good  choice  because  it  is  transparent  in  the  visible  and  IR  regions  of  the 
spectrum  (through  7  microns)  [1 1 :  38, 5:  21].  The  only  drawback  of  using  this  material  is 
its  low  thermal  diffusion  coefficient,  which  puts  a  restriction  on  how  fast  the  IR  beam  can 
be  modulated.  If  the  beam  receives  light  at  too  high  a  frequency,  no  equilibrium  will  be 
reached.  Measurement  of  deflection  is  made  during  steady  state  operation,  where  the 
buildup  of  heat  in  the  top  layer  equals  the  amount  leaving.  No  one  has  yet 
mathematically  characterized  the  mechanical  behavior  of  multi-layer  reeds.  Dr  William 


2.2 


P.  Baker,  Associate  Professor  of  Mathematics  at  the  Air  Force  Institute  of  Technology 
has  recently  developed  a  description  of  the  thermal  behavior  of  the  three-layer  reed. 

Static  Mechanical  Analysis 

The  fundamental  quantities  used  in  describing  the  mechanical  motion  of  the  reed 
are  given  below.  Additional  parameters  of  interest  will  be  defined  as  they  are  discussed. 
Only  the  silicon  and  aluminum  layers,  layers  two  and  three  respectively,  are  considered 
in  the  static  analysis  of  the  reed. 


Table  2.1:  Definition  of  Terms 


Name 

Ouantitv 

Name 

Ouantitv 

z 

Vertical  deflection  (m) 

m 

Sensor  Temperature  (K) 

s 

Sensitivity  (m/W) 

K: 

Young’s  Modulus  (unitless) 

«2,  a3 

Thermal  Expansion  Coefficients  (1/K) 

m 

Length  of  the  cantilever  reed  (m) 

^2,  X.3 

W: 

Width  of  the  cantilever  reed  (m) 

t2,t3 

Thickness  (m) 

mm 

Power  absorbed  by  reed  (W) 

tm 

Ambient  Temperature  (K) 

Numerical  values  for  each  material  can  be  found  in  the  Appendix  A.  Ambient 
temperature  is  that  temperature  at  which  the  beam  lays  flat. 

Current  literature  discusses  two-layer  cantilevers  and  gives  the  theory  for  the 
static  case,  where  the  reed  has  previously  absorbed  some  amount  of  energy  and  has 
deflected.  Assuming  that  heat  losses  due  to  convection  and  radiation  are  negligible,  the 
reed  sensor  can  be  approximated  by  a  two-layer  cantilever  acting  under  uniform 
(throughout  the  volume)  heating,  with  an  aspect  ratio  (length  /  width)  of  at  least  ten  and  a 
length  24  times  greater  than  thickness  [15:  83].  This  last  condition  is  always  met.  If  the 
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required  aspect  ratio  is  met,  deflection  occurs  only  as  a  function  of  x  as  defined  in  Figure 
2.1.  Deflection  of  reeds  with  a  smaller  aspect  ratio  depends  also  on  the  width  axis.  The 
multi-mode  mechanical  analysis  necessary  to  treat  these  cases  is  rather  complicated  and 
has  only  been  developed  for  deflections  of  single  layer  cantilevers  [12:  2480]. 

Thickness  of  the  silicon  layer  is  not  variable,  as  these  are  prefabricated,  but  the 
aluminum  thickness  can  be  varied.  The  sol-gel  thickness  could  also  be  varied  using 
numerical  calculations,  but  was  held  constant  for  this  modeling.  Current  reeds  are  made 
with  a  thick  silicon  layer  and  a  thinner  aluminum  layer.  As  the  thickness  of  the 
aluminum  layer  is  increased,  deflection  of  the  reed  sensor  will  also  increase  [4:  32].  The 
sensitivity  peaks  at  some  optimum  value.  The  sensitivity  as  a  function  of  the  ratio  of 
aluminum  thickness  to  silicon  thickness,  n,  is  [6:  1 14]: 

2a2L3  Q3-l)(tt  +  l)  (2-1) 

A ^2  W  K(n)(yn  +  \) 

where 

K{n)  =  4+6ra  +  4n2  +0w3  +— . 

<j>n 


Here  <|)  is  the  ratio  of  the  elastic  modulus  of  aluminum  to  that  of  silicon,  y  is  the 
ratio  of  thermal,  and  P  is  the  ratio  of  thermal  expansion  coefficients.  This  equation  can 
be  used  to  determine  the  optimum  thickness  of  the  aluminum  layer  in  terms  of  the  silicon 
thickness. 

By  modeling  the  reed  as  a  lumped-mass  spring  system,  a  restorative  force  (spring) 
constant,  kc,  and  an  effective  mass  of  the  reed,  meff,  can  be  used  to  calculate  mechanical 
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resonant  frequency.  The  following  mechanical  resonance  equations  are  appropriate  for 
rectangular  bi -layer  beams  with  an  aspect  ratio  of  at  least  8  [13:  13,  23]: 


kc  —  ~{E3t\  +  E2t\) 

4  L 

meff  =0.24  wL(p3t3+p2t2) 


co 


O 


t2  |  (E3n2  +  E2) 
2L2\24(p3n  +  p2) 


(2.2) 


For  each  layer,  i,  density  is  represented  by  pi,  elastic  modulus  by  Ej,  and  thickness 
by  tj.  Note  that  with  increased  mass  of  material,  the  resonant  frequency  will  go  down. 

Solution  to  the  Thermal  Model 

During  this  research,  Dr.  Baker  developed  equations  describing  heat  transfer  in 
the  reed’s  layers.  They  are  based  on  the  fundamental  heat  transfer  equation,  written  for 
each  layer  of  the  reed.  A  complete  derivation  can  be  found  in  the  Appendix  B.  The 
origin  is  placed  at  the  far-left  comer  of  the  reed,  between  the  first  and  second  layers,  as 
seen  below  in  Figure  2.1. 


y 


hi 

com _ 

Sol-gel 

h2 

Silicon 

7 

h3] 

_ 

Aluminum 

Figure  2.1:  Schematic  of  the  Reed 
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Heat  transfer  equations,  written  for  all  three  layers,  must  be  solved  to  determine 
the  time  based  temperature  profile.  Here  vertical,  y,  and  horizontal,  x,  directions  are 
modeled.  The  assumption  is  made  that  the  temperature  profile  does  not  vary  with  reed 
width.  Convection  and  radiation  remove  a  negligible  amount  of  heat  from  the  reed  and 
are  ignored.  Energy  is  deposited  in  the  first  layer  by  a  modulated  IR  beam  at  a  constant 
intensity  for  half  of  a  modulation  period.  This  development  assumes  uniform  energy 
deposition  in  the  volume  of  first  layer.  Also  assumed  is  that  the  temperature  at  the  left 
end  of  the  reed  is  held  constant  at  ambient  temperature  by  the  clamp  mechanism.  This 
temperature  is  also  the  initial  temperature  of  the  reed.  The  other  three  boundaries  of  the 
reed  are  taken  as  insulated. 

The  solutions  for  each  layer  are  coupled  by  the  additional  boundary  conditions  of 
continuity  of  temperature  and  flux  across  the  interfaces.  Let  w(4,r|,x)  represent  a 
temperature  fluctuation  from  the  initial,  and  ambient,  temperature.  Here  the  spatial 
dimensions  have  been  scaled  to  the  overall  length  of  the  reed,  that  is  length  x  =  L%,  and 
depth  y  =  J]L.  The  thickness  of  each  individual  layer  is  also  scaled  with  respect  to  the 
overall  length  resulting  in  hx  =  Loti.  Time,  t,  is  scaled  to  the  period  (T)  of  the  modulated 
IR  light  beam,  or  t  =  Tt.  The  heat  transfer  equation,  describing  the  fluctuation  of 
temperature  from  the  ambient  and  written  for  each  layer  i,  is: 


=  diV1w{l)  (£,r\,T)+8iXF(r)  . 


Here,  energy  deposition  rate  from  the  IR  beam  is  represented  as  F(x)  and  is  present  only 
in  the  first  layer,  hence  use  of  the  Kroneker  delta,  8j,i.  The  specific  function  of  this 
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energy  source  will  be  considered  later.  The  term  dt  in  equation  (1)  is  a  scaled  diffusion 
term  which  includes  some  constants  as  a  result  of  scaling: 


Di  is  the  thermal  diffusivity  for  each  layer  i.  The  left  end  of  the  reed  is  held  at 
constant  temperature,  so  the  temperature  fluctuation,  w,  is  always  zero  at  %  =  0.  The 
reed’s  top  and  bottom  are  insulated,  meaning  there  will  be  no  temperature  variation  in  the 
vertical  (y)  direction  at  the  top  of  the  sol  gel,  T|  =  oq,  nor  at  the  bottom  of  the  aluminum 
layer,  r)  =  -(a2  +a3).  There  is  no  temperature  variation  along  the  horizontal  (x)  at  the 
insulated  right  end  of  the  reed,  £  =  1. 

Again,  the  separate  equations  are  coupled  by  continuity  of  temperature  across  the 
boundaries  and  continuity  of  heat  flux  across  the  boundaries.  Finally,  the  reed  starts,  in 
an  undeflected  position,  at  the  ambient,  initial  temperature:  w(i)(£,t7,0)  =  0 . 

To  satisfy  boundary  conditions  at  §=  0,  1,  the  solution  will  have  the  form: 

w(l)  (£,r?,T)  =  ^qk(i)(ri,r)Sin(AkZ)  y~A) 

k=  0 

where 

Ak  =  (2k +1)—  for  k  =  0,1,2.. 

<2 

and  qkfl)  are  mode  shapes  to  be  determined. 

Substituting  this  solution  into  the  original  differential  equation,  and  Laplace 

transforming  the  result,  yields  an  equation  in  Qk\ri,s)  =  L|^0  (r/,t)j.  The  solution  is 

easily  obtained.  Directly  inverting  the  Laplace  transform,  via  calculus  of  residues, 
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requires  the  roots  in  the  transform  domain  of  an  algebraic  equation  (see  Appendix  B). 
These  simple  roots,  labeled  (-ak ),  lie  along  the  negative  real  axis,  and  represent  the  rate 
of  diffusion  through  the  reed.  The  convolution  theorem  is  used  to  determine  the 
q(,\ri,r) ,  once  the  incoming  energy  is  determined.  The  thermal  deposition  rate  from  the 
IR  source,  F(t),  is  modulated  by  a  chopper,  and  takes  the  form  of: 

F(r)  =  kaE0(r)  =  ka 


Here  t0  represents  the  fraction  of  a  period,  T,  that  the  beam  is  turned  on.  Only  a  portion, 
ka,  of  the  incident  energy,  Eo,  is  actually  converted  to  heat  present  in  the  sol-gel  layer. 
The  parameters  ka  (absorption  of  the  sample  trapped  within  the  sol-gel  layer)  and  I0, 
(intensity  of  the  incoming  energy  beam)  can  be  determined  experimentally. 

Now  Ik(t)  is  defined  by  Ik  ( t )  =  {e~Gk^~T^F(t)dz  and  is  found  to  be 

J  0 


40 0 


kq[o  e~ffkr  (e  11  (eakP  1)  eakP  |  KJo_  T 

°k  (e°k  -i)  J  ak  \e-ok(.(r-p)-t0)  t 

0  rr, 


(2.6) 


The  first  term  in  /*( t)  in  Equation  2.6  accounts  for  the  build  up  and  decay  of 
energy  for  each  of  the  previous  cycles.  The  other  accounts  for  the  cycling  of  the  current 
period,  p.  Heat  in  the  top  layer  builds  up  while  the  beam  is  on.  Once  the  beam  has  shut 


off  (t/T  >  to),  the  temperature  slowly  decays,  as  heat  is  transferred  to  the  other  layers  and 
the  clamp. 
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Recall  that  (-ok )  are  roots  of  the  solution  in  the  transform  domain,  needed  to 


invert  the  Laplace  transform. 


Phase  analysis  can  be  obtained  from  the  model  by  noting  the  time  it  takes  for  the 
temperature  to  rise  and  fall  to  a  defined  point.  A  convenient  point  is  half  of  the 
magnitude  of  temperature  change. 


Summary  of  the  Thermal  Solution. 

The  final  solution  for  each  layer  is  comprised  of  a  series  of  solutions.  Each  series 
term,  k  includes  a  dynamic  description  of  the  energy  deposition  (Equation  2.6), 

modulated  in  the  horizontal  (x)  direction  by  </>k(£)  =—Sin(Xk%)  . 

K 

The  temperature  diminishes  as  it  passes  through  the  depth  (y)  of  the  reed.  This 
behavior,  of  the  scaled  temperature,  is  shown,  for  each  layer  i,  by  y/‘k(ri): 


i 

YM  Cos(  -j==  (tj  of j )) 
V“i 


(i) 


(2.7) 


»(2) 


(2) 


9  0), 

Yk  (V)  =  AkCosh(—j=  (1 rj +a2))  +  Bk  Sinh(-f=  (r\  +  a2 )) 

V“  2  V^2 


where 


Ak  =Cosh((Qik)Yi) 


B„=  &L 


d ,  (o 


(3) 


d2  o) 


(2) 

k 


Sinh(a)(k3)y3) 


(0 


(3) 


Yk(rl)  =  Cosh(-~=  ( rj+a2+  a3 )) . 


The  second  layer  is  far  more  complex,  being  trapped  between  a  layer  that  is 
constantly  introducing  more  energy  into  the  system  (the  top)  and  a  layer  that  quickly 
equilibrates  any  temperature  differences  it  sees. 

For  each  unit  of  energy  deposited,  there  is  a  temperature  change: 

c f  =  7  (1K2~  \eoPsinh(co^Y2)Cosh(co^y3)  +  ^ojf)Cosh(cof)Y2 )Sinh(o)f>y3 )  (2'8) 

(ffl*  )  Dk  l  Vrf2 

c(2)  _  M3)  _  2  pT Sinjaln) 
k  k  id2  co^Dk  ' 

The  function  Dk  is  defined  in  Appendix  B.  For  convenience,  the  following  is 
written: 


^zk+dxX\  =  io)(k} 

Yi  and  ^zk  +  d2Xk  =cok^ . 

~Jzk+  d-itfk  = 

The  solution  to  the  thermal  behavior  for  a  three-layer  reed  is  then: 


=  ±cflk  (t  )<pk  (^)^°  (??) 


k=0 


(2.9) 


Dynamic  Mechanical  Modeling 

The  differential  equation  for  the  two-layer  beam  is  [15:  1 18, 2:  3793]: 

7  2  „ 


4# =  6(a2  -ax)hll2{T.T0){x,t) 


(2.10) 


dx 


t\K 


The  reed  is  oriented  such  that  the  end  at  x  =  0  is  fixed  and  the  other,  at  x  =  L  is 
free.  The  fixed  end  shows  no  vertical  movement  and  no  slope.  Replacing  the  term  (T- 
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T0)(x,t)  with  Baker’s  thermal  solution,  and  solving  for  the  deflection  at  the  top  of  the 
reed’s  free  end  results  in, 


zx=l 00  =  6(«i  -  a2 )  tXY3u  E  cfh  (J)fk  0-)vil)  («i ) 


(2.11) 


tZK 


k=Q 


where  the  functions  in  the  series  are  defined  as  before,  with  the  exception  of: 

At  A/, 


(2.12) 


As  mentioned  earlier,  there  are  roots  associated  with  each  term  of  the  series. 
These  roots  can  be  viewed  as  an  effective  thermal  diffusion  constant  for  each  mode.  The 
roots  for  each  series  of  equations  must  be  re-calculated  for  each  parameter  change.  This 
makes  it  difficult  to  optimize  for  length,  thickness,  and  operating  frequency.  To 
circumvent  this  problem,  and  obtain  deflection  as  a  function  of  length  and  chopper 
frequency,  several  approximations  were  required. 

Once  the  deflection  is  known,  the  signal-to-noise  ratio  (SNR)  can  be  determined, 
as  it  is  dependent  on  deflection.  This  analysis  assumes  thermal  noise  dominates  all  other 
sources.  A  two-layer  reed  sensor,  modeled  as  a  spring  mass  system,  will  have  resonant 
frequency  <a„  and  spring  constant  kc,  when  operating  at  ambient  temperature  T.  It  will 
have  quality  factor  Q,  which  can  be  determined  by  the  shape  of  the  first  resonance  peak. 
Finally,  the  signal  to  noise  ratio,  SNR,  is  [13:  31]: 


SNR  —  zx=i . 


\Q%Oo 

4  kbT 


(2.13) 


where  kb  is  the  Boltzmann  constant. 

Maximizing  the  SNR  is  the  basis  for  optimizing  the  reed  sensor’s  geometry. 
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3.  Methodology 


Experimental  Procedure. 

Preparation  of  Reed  Sensor.  Commercially  prepared  silicon,  with  a  thickness  range  of 
five  to  ten  microns,  is  cut  to  size,  and  an  evaporator  is  used  to  coat  a  thin  layer  of  gold  on 
one  side.  A  liquid  solution  containing  the  sol-gel  materials  is  dropped  onto  the  gold  layer 
and  “spun”  in  a  spinner  to  achieve  a  sub-micron  thick  film.  Drying  the  new  reed  removes 
any  moisture  and  some  organic  impurities  present  in  the  sol-gel  solution  [8:  259]. 

Finally,  the  aluminum  is  evaporated  onto  the  other  side  of  the  silicon.  The  process  takes 
about  10  days,  and  has  been  nearly  perfected  by  AFIT  technician  Ms.  Belinda  Johnson 
and  Dr.  Guangming  Li. 

Once  ready  for  use,  a  reed  is  clamped  on  a  silicon  base  and  inserted  into  the  AFM. 
The  base  is  represented  in  Figure  1 . 1  and  consists  of  one  layer  secured  in  the  AFM,  and 
two  small  pieces  of  silicon  at  either  end.  The  first  piece  acts  as  a  spacer,  keeping  the  reed 
from  touching  the  base.  The  second,  as  a  stopper  at  the  end  of  the  reed,  helps  damp 
unwanted  vibrations.  As  a  bonus,  this  became  a  suitable  platform  to  measure  reed 
deflection,  as  the  force  of  the  AFM  tip  will  continue  to  push  the  reed  down,  if  there  is 
nothing  stable  below.  The  AFM  can  then  measure  the  upward  force  caused  by  a  change 
in  temperature  of  the  reed  once  it  has  a  stable  platform  against  which  to  work. 

Testing  with  the  AFM.  The  AFM  was  used  to  measure  the  deflection  response  to 
changes  in  chopper  frequency,  wavelength  of  light,  sample  concentration  and  use  of  a 
reflector.  Unfortunately,  small  variations  in  the  starting  position  of  the  AFM  tip  prevent 
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the  use  of  software  to  subtract  out  such  background  and  zero  signals.  A  schematic  of  the 
laboratory  setup  is  given  below: 


Figure  3.1  Laboratory  Setup 

There  are  several  ways  to  measure  phase  and  amplitude  information.  A  lock-in 
amplifier,  the  standard,  but  slow,  measurement  device  in  this  field,  provides  accurate 
phase  information,  but  amplitude  measurements  are  highly  noisy,  due  to  the  laboratory 
environment.  A  much  faster  method  to  measure  amplitude  is  found  in  the  software 
accompanying  the  AFM.  Image  processing  software  is  able  to  correct  for  any  tilt  of  the 
reed,  perhaps  owing  to  variations  of  the  base  and  stopper  heights.  The  software  can  also 
approximate  a  low  pass  filter,  smoothing  out  the  deflection  data  and  removing  power 
spikes  that  often  threw  off  the  lock-in  measurements.  Phase  information  cannot  be 
determined  using  the  AFM  without  reference  to  the  chopper  frequency  signal. 
Measurement.  The  light  used  in  most  of  this  analysis  was  in  the  visible  range.  The 
gratings  used  were  blazed  at  500  nm  and  were  relatively  inefficient  in  the  300-400  nm 
range,  but  efficient  from  500-650.  The  absorber  used  in  this  study  was  a  laser  dye, 
pyridine,  which  had  two  absorption  peaks:  A  major  peak  at  485  nm  and  a  much  smaller 
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peak  at  about  350  nm.  The  possibility  of  attributing  deflection  peaking  from  the  light 
source  to  that  resulting  from  the  sample  is  reduced  by  taking  source  spectra  to  determine 
the  amount  of  light  absorbed  by  the  sample  for  each  wavelength. 

Determination  of  Parameters.  The  mathematical  model  will  show  the  trends  in  reed 
behavior,  but  will  not  provide  specific  energy,  deflection,  or  absorption  predictions 
without  calculating  ka,  the  absorption  of  the  sample  in  the  sol-gel  layer  and  Io,  the 
intensity  of  the  incoming  energy  beam.  An  optical  power  meter  is  used  to  determine  the 
fraction  of  original  energy  from  the  light  (90  to  100  Watts)  that  actually  hits  the  reed.  The 
next  step,  measuring  Io,  requires  measuring  the  absorption. 

Mathematical  Modeling. 

Software  Package.  Mathematical  modeling  was  done  using  the  Mathematica  software 
package.  Mathematica  was  chosen  mainly  for  its  graphic  capabilities.  Compared  to 
other  popular  software  packages,  it  is  more  user  friendly  than  Fortran,  and  is  able  to 
calculate  multiple  series  and  present  three-dimensional  graphics  with  far  less 
manipulation  than  MathCad.  Three-dimensional  graphics  were  not  heavily  used  in  this 
research,  but  it  is  expected  that  someone  will  take  over  the  research  and  derive 
relationships  involving  thickness  of  the  top  layer,  length,  and  frequency.  These  new 
relationships  will  produce  parameter  values  for  an  optimum  reed,  and  the  relationships 
will  be  easier  to  see  with  three-dimensional  snapshots.  This  software  package  does  have 
a  distinct  disadvantage  when  working  with  units,  as  can  be  seen  in  the  code  presented  in 
Appendix  C. 
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Approximations.  Applying  the  three-layer  thermal  solution  to  the  two-layer  mechanical 
solution,  required  use  of  several  approximations  and  assumptions.  First,  because  the 
expansion  of  the  second  and  third  layers  determines  the  deflection  of  the  reed,  the 
temperature  of  the  second  layer  along  its  top  was  used  to  describe  the  temperature  profile 
of  the  reed.  This  profile  fed  into  the  deflection  equation  as  presented  in  the  previous 
chapter.  This  approach  was  later  validated  when  the  bottom  two  layers  were  found  to 
exhibit  the  same  temperature  profile. 

Parameter  space  exploration  was  made  difficult  by  the  presence  of  the  root  terms, 
which  needed  to  be  re-calculated  for  each  specific  reed  and  chopper  frequency.  To  make 
the  optimization  more  tractable,  approximations  were  found  for  the  root  functions  (-crk ) 

and  the  derivative  functions  ( Dk ).  Approximations  for  both  the  roots  and  the  derivative 

terms  were  made  relative  to  the  baseline  case  of  a  reed  0.01  meters  in  length  and  operated 
at  1  Hz  (dc).  These  functions,  introduced  in  the  derivation  of  the  thermal  solution,  were 
found  to  follow  the  simple  ratio  relationships  presented  in  the  fourth  chapter. 

Inputs.  The  Mathematica  code  produced  during  this  project  is  presented  in  Appendix  C. 

It  requires  user  knowledge  of  the  reed  parameters  such  as  length,  thickness  and  chopper 
frequency,  as  well  as  material  properties  such  as  density,  elastic  modulus,  thermal 
conductivity,  thermal  expansion,  and  heat  capacity  for  each  layer.  Further  required 
quantities,  diffusivity  for  example,  can  be  calculated  from  those  given  above. 

Outputs.  The  code  calculates  first  the  thermal  profile  of  a  given  reed  and  the  uses  that 
profile  to  determine  a  dynamic  deflection  profile  for  the  reed  under  consideration.  The 
static  portion  of  the  model  can  be  used  to  explore  resonant  frequency,  optimum  thickness 
and  ratio  of  the  bottom  two  layers  of  the  reed. 
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Typical  Values.  The  code  has  been  adapted  to  be  useful  over  a  wide  range  of  reed 
lengths  and  chopper  frequencies.  Most  analysis  was  done  for  the  same  case:  a  reed  2.25 
cm  long,  with  the  thickness  of  the  three  layers  (in  descending  order)  being  1  micron,  8 
microns,  0.5  micron.  The  chopper  frequency  was  typically  set  at  8  hertz.  These 
parameters  were  used  because  they  were  the  parameters  of  the  typical  reeds  used  in  the 
laboratory.  The  absorption  coefficient,  ka,  was  not  determined,  and  left  as  unity  in  many 
calculations.  The  degradation  of  the  light  signal  was  accounted  for  by  using  values 
measured  in  the  lab  [7:  1]. 
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4.  Results 


Model  Results 

Optimum  Thickness  Ratio.  In  Chapter  2  we  discussed  the  variation  in  sensitivity  as  a 
function  of  n,  the  ratio  of  silicon  and  aluminum  thickness.  Figure  4.1  below  is  a  plot  of 
sensitivity  as  a  function  of  the  ratio  n  and  reveals  an  optimum  ratio  of  n  =  0.325.  Thus, 
for  a  silicon  thickness  of  8  microns,  the  optimum  aluminum  thickness  is  2.6  microns. 

The  aluminum  layer  cannot  be  made  that  thick,  as  current  processing  techniques  limit  the 
thickness  to  half  a  micron.  Beyond  this  thickness,  because  of  the  curling  following  A1 
evaporation  at  elevated  temperatures,  the  reed  assumes  a  curled,  rather  than  flat  shape  as 
desired. 
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Figure  4.1  Sensitivity  as  a  Function  of  (n)  From  Equation  2.1 
Showing  Optimum  Thickness  Ratio 

The  silicon  thickness  is  not  as  constant  as  desired-  being  found  in  the  lab  to  range 
anywhere  from  four  to  sixteen  microns.  The  optimum  value  for  aluminum  remains  0.325 
regardless  of  the  silicon  thickness,  another  result  of  static  modeling.  One  way  to  increase 
performance  of  the  reed  sensor  is  to  use  methods  that  will  allow  for  a  larger  aluminum 
thickness. 
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Resonant  Frequency.  Static  two  layer  mechanical  modeling  is  used  to  calculate  resonant 
frequency  for  a  reed  of  silicon  and  aluminum.  The  silicon  layer  is  assumed  8  microns 
thick,  and  the  aluminum  layer  takes  on  values  from  0.5  to  8  microns  thick.  The  ratio  of 
aluminum  to  silicon  thickness  is  noted  as  n,  according  to  the  resonance  equation, 

Equation  2.2.  The  resonant  frequency  has  a  minimum  at  where  the  aluminum  thickness 
=  0.6  of  the  silicon  thickness.  This  is  far  too  large  given  current  reed  production 
procedures.  Other  methods  of  bringing  the  resonance  frequency  down  to  reasonable 
chopper  frequencies  must  be  explored. 

The  Thermal  Solution 

Predictions  obtained  for  the  dynamic  thermal  behavior  of  the  reed  sensor  met 
expectations.  There  is  a  relatively  large  temperature  gradient  in  the  first  layer,  and  the 
bottom  two  layers  show  almost  constant  temperature.  At  the  end  of  the  reed,  the 
temperature  profile  is  constant  throughout  all  three  layers.  This  condition  actually  occurs 
long  before  the  end  of  the  reed,  at  a  point  slightly  less  than  half  way  down.  The 
temperature  increases  nearly  linearly  but  soon  reaches  a  steady  state  condition  where 
each  layer’s  temperature  oscillates  between  two  constant  values.  It  is  this  minor  vibration 
that  is  measured  as  deflection  due  to  light  absorption.  This  steady  state  is  reached  in  less 
than  100  periods.  In  figure  4.2  the  rise  to  equilibrium  can  be  seen.  In  figure  4.3,  after 
200  cycles  of  the  light  beam,  we  see  the  steady  state  operation  of  the  reed  oscillation. 
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Figure  4.2:  Thermal  Response  of  Multi-Layer  Reed  in  Time 
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Figure  4.3:  Steady  State  Thermal  Behavior  of  Multi-Layer  Reed 


During  steady  state  operation,  the  temperature  of  the  reed  just  before  the  start  of 
any  period  (just  before  the  beam  turns  on)  is  constant.  The  difference  in  temperatures  can 
be  seen  in  Figure  4.4,  which  is  a  plot  of  the  temperature  fluctuation  as  a  function  of  depth 
into  the  reed.  This  snapshot  is  taken  halfway  through  the  2001st  period,  just  as  the  beam 
shuts  off.  At  the  beginning  of  the  period,  the  temperature  profile  had  the  same  shape,  but 
the  variation  was  so  slight  that  it  would  appear  as  a  straight  line,  if  plotted.  The 
temperature  profile  is  nearly  constant  when  the  light  begins  to  strike  the  beam,  but  a 
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difference  in  profile  throughout  the  depth  of  the  reed  quickly  develops,  peaking  when  the 
light  turns  back  off.  The  calculation  was  done  near  the  clamped  end  of  the  reed  to  show 
the  maximum  temperature  difference  between  the  layers. 


Figure  4.4:  Temperature  Change  vs.  Depth  as  Beam  Turns  Off, 

Steady  State  Operation 

The  difference  in  temperature  between  the  layers  is  about  10'10  Kelvin  for  the 
baseline  case  mentioned  in  Chapter  3.  The  difference  in  temperature  between  periods  is 
still  small,  on  the  order  of  10'8  Kelvin.  This  seemingly  small  temperature  change 
corresponds  to  a  calculated  .01  nm  deflection  for  the  typical  four-layer  reed.  Measured 
deflection  values  for  this  reed  ranged  between  30-200  nm. 

Thermal  Solution  Applied  to  Mechanical  Model.  Recall  that  the  thermal  solution  was  an 
infinite  series,  each  term  of  which  containing  a  root  which  described  the  overall  heat 
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transfer  of  the  reed  for  that  mode.  These  roots  had  to  be  numerically  recalculated  with 
each  change  in  reed  geometry.  An  analytic  representation  of  these  quantities  would  allow 
for  more  convenient  application  of  the  thermal  and  mechanical  solutions  to  the  reed 
sensor  problem. 

To  find  such  a  relationship,  these  roots  were  calculated  for  several  different  reeds, 
all  compared  to  a  baseline  reed  0.01  meters  in  length  modulated  at  1  Hertz..  First, 
plotting  the  root  corresponding  to  each  series  term,  k,  revealed  a  parabolic  relationship  in 
k.  The  length  was  then  varied,  and  the  roots  were  observed  to  be  proportional  to  one  over 
square  of  length  of  the  reed.  Finally,  chopper  frequency  was  varied  and  the  roots  varied 
as  one-over-frequency. 

These  three  effects  were  combined  by  first  determining  the  parabolic  constants  for 
the  roots  as  a  function  of  k  for  the  baseline  case.  That  root  value  is  then  multiplied  by  the 
appropriate  ratios  in  L  and  to  to  determine  the  value  of  the  root  for  other  lengths  and 
frequencies.  The  roots  (-  <7k )  are  then  represented  by: 

ak  ( k,L ,  co)  =  (-2.02072  +  8.1 9207  *  k  -  8.202 16  *  k 2 ) . 

coif 

The  term  .0001  has  the  units  Hertz-Square  meters,  and  so  the  scaled  nature  of  the 
problem  can  be  maintained  by  inputting  L  in  meters  and  to  in  Hertz.  This  term  stems 
from  proportions  of  L  and  CD,  that  is,  the  values  L=  .01  (roots  vary  inversely  as  L  squared) 
to  =  1  Hz  used  in  calculating  the  original  parabolic  constants. 

The  derivative  terms  ( Dk ),  introduced  in  Chapter  2  were  also  approximated. 

They  too  were  found  to  obey  simple  relationships,  but  did  not  show  significant  variance 
with  L.  Dk  ( k,(0)  =  (.00099027 15— 3 .32884 1  *  10-9 -fe — 7.0945763 16*1 0 “10  £2 )  V© 
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The  k  and  fc2  constants  above  may  seem  relatively  small,  but  keep  in  mind  that  the 
derivative  terms  are  just  as  small  (on  the  order  of  10E-4).  These  small  constants  make  a 
large  difference  in  approximating  higher  series  terms.  The  thermal  code  itself  does  not 
require  use  of  series  terms  higher  than  100,  but  the  after  approximating  the  roots  and 
derivatives  in  this  manner,  it  may  became  necessary  to  run  the  series  out  to  about  200, 
especially  for  smaller  reeds.  Ignoring  these  tiny  terms  can  result  in  a  30%  error  in  root 
terms. 

Table  4. 1  shows  the  error  introduced  when  using  the  data  fits.  The  Data  Point 
column  shows  the  conditions  under  which  the  root  was  calculated  as  follows: 

Root  [Length  in  Meters,  Chopper  Frequency,  Series  Term]. 

For  example,  the  first  entry  represents  the  root  calculation  for  a  reed  that  was  .03 
meters  long,  the  chopper  frequency  was  8  Hertz,  and  the  thermal  series  was  taken  out  12 
terms.  A  similar  table  highlighting  the  error  in  the  derivative  terms  is  in  Appendix  A. 


Table  4.1:  Error  in  Roots  with  Approximation 


Data  Point 

Actual  Root 

Approximated  Root 

%Error 

Root[.03, 14,12] 

-10.1722 

-10.1721 

.001% 

-131.696 

Root[.01, 1,500] 

-2.007825* 10A6 

-2.004645* 10A6 

1.9% 

Root[.004,29,250] 

-106512. 

-110041. 

3.3% 

The  variance  of  the  roots  and  derivatives  with  respect  to  length  and  frequency  are 
based  on  the  change  from  an  initial  value,  calculated  for  a  reed  one  centimeter  in  length 
modulated  at  1  Hz.  These  approximations  allow  thermal  profile  to  be  determined  as  a 
function  of  length  and  chopper  frequency,  and  have  small  error  compared  to  the 
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predictions  obtained  by  Baker’s  original  equation.  The  error  does  grow,  however,  for 
very  small  reeds  (less  than  half  a  centimeter)  at  very  high  series  term  numbers.  The  roots 
at  this  point  are  so  large  that  this  error  hardly  effects  the  overall  temperature  prediction. 

Overall  predictions  in  temperature  showed  no  greater  than  1%  difference.  An 
important  quantity  is  the  difference  in  temperature  that  the  reed  experiences  during  steady 
state  operation.  The  approximations  modeled  this  temperature  change  with  an  error  less 
than  .01%.  This  study  varied  the  reed  length  from  four  millimeters  to  25  centimeters  and 
ran  the  chopper  frequency  from  one  Hertz  up  to  about  500  Hertz. 

The  model  calculated  the  expected  temperature  profile  across  the  length  of  the 
reed  (x  direction ),  as  shown  below  in  Figure  4.5.  There  are  three  curves  shown  in  each 
graph:  The  lowest  curve  shows  the  reed  after  only  one  period’s  worth  of  light  has  been 
absorbed.  The  middle  shows  the  temperature  profile  after  ten  periods  have  passed,  and 
the  highest  curve  shows  the  steady  state  condition  resulting  after  100  periods  have 
passed.  This  is  temperature  in  the  top  layer.  The  deflection  of  the  reed  will  be  seen  in 
time  variations  of  the  height  of  that  steady  state  curve.  The  thermal  solution  and  the 
approximated  solution  produce  the  same  curves. 
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Figure  4.5:  Temperature  Profile  Along  Reed  Length 


Phase  Lag 

Phase  lag  was  explored  by  noting  comparing  the  expanse  of  time  between  rise  and  fall  to 
a  specific  temperature  point.  The  time  at  which  the  reed  first  heats  up  to  this  specific 
temperature  is  labeled  t-early,  4,  and  the  time  to  fall  back  is  labeled  t-late,  tj.  The  change 
in  the  difference  between  these  two  times  for  each  series  term  quickly  drops  to  a  constant 
value  of  one  half  (See  Appendix  C).  This  suggests  that  accurate  phase  information  can 
be  determined  with  only  the  first  series  (mode)  term. 

Below  is  a  plot  showing  combining  thermal  output  of  the  first  layer  from  two 
different  times.  The  first  curve  represents  the  dynamic  temperature  profile  during  the 
first  period.  The  second  was  calculated  during  steady  state  conditions.  Phase  in  the 
signal  is  seen  in  the  difference  in  the  rise  and  fall  of  the  curve. 
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Signal  to  Noise  Ratio  (SNR').  Mathematical  modeling  produced  expected  results  when 
tracking  the  effect  of  length  and  frequency  on  the  overall  deflection  of  the  reed.  As 
length  was  increased,  total  deflection  of  the  reed  and  noise  present  increased.  As  the 
frequency  was  increased  (beyond  mechanical  resonance  frequency  of  the  reed)  the 
deflection  went  down. 

The  model  (presented  in  Appendix  C)  also  shows  expected  results  in  that  the 
longer  the  reed  gets,  the  less  sensitive  it  becomes.  Although  the  deflection  of  a  reed 
increases  along  with  its  length,  the  longer  reeds  have  very  low  resonance  frequencies. 
These  low  frequencies  result  in  a  reed  sensor  that  is  very  sensitive  to  low  frequency 
background  noise.  SNR  was  predicted  to  peak  at  the  resonant  frequency  and  then  drop  off 
sharply.  The  resonant  frequency  of  the  long  reeds  (about  2  cm  in  length)  was  as  low  as 
our  lowest  operating  frequencies  (8Hz)  and  so  resonance  peaking  was  not  observed.  A 
three-dimensional  plot  of  predicted  SNR,  as  well  as  the  calculations  that  produced  it,  is 
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presented  in  Appendix  C.  For  each  length,  the  maximum  SNR  will  be  obtained  by 
operating  just  short  of  the  mechanical  resonance  frequency  of  the  reed. 

Experimental  Results 

Long  Vs  Short  Reed  Sensors.  As  expected,  there  was  a  distinct  difference  in  the  signals 
obtained  from  small  reeds  and  long  reeds.  The  long  reeds  (two  centimeters  and  longer) 
were  far  more  susceptible  to  acoustic  noise,  while  the  short  reeds  (1.5  centimeters  and 
longer)  produced  a  zero  signal  except  at  very  low  frequencies  (due  to  low  frequency 
noise)  and  near  the  resonance  frequency.  The  resonant  frequency  of  the  long  reeds 
(about  2  cm  in  length)  was  as  low  as  our  lowest  operating  frequencies  (8Hz)  and  so 
resonance  peaking  was  not  observed.  Resonance  was  observed  for  short  reeds,  as  will  be 
shown  later. 

Resonance  Frequency.  The  resonant  frequency  calculations  assume  that  the  reed  meets 
beam  theory  requirements  presented  in  Chapter  2.  To  see  the  effect  on  reeds  which 
violate  these  requirements,  a  reed  1.25  centimeters  in  length  and  3  millimeters  wide  was 
prepared  and  tested  for  its  deflection  response  to  chopper  frequency.  Because  the  aspect 
ratio  is  far  from  ten  (being  4.5),  some  deviation  from  the  predicted  resonance  of  330 
Hertz  was  expected.  Not  only  did  the  resonance  peak  appear  at  about  425  Hertz,  but  also 
the  shape  of  the  deflection  curve  revealed  higher  order  resonances  present.  This  data  is 
presented  as  below  as  Figure  4.7.  Arrows  are  placed  over  the  expected  and  actual 
resonance  peaks.  The  aspect  ratio  is  so  low  that  beam  mechanics  no  longer  apply.  The 
deflection  is  no  longer  one-dimensional.  There  is  now  resonance  behavior  across  the 
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width  of  the  reed.  To  account  for  this  kind  of  vibration,  the  thermal  and  mechanical 
solutions  must  be  derived  using  plate  mechanics. 


Chopper  Frequency  (Hertz) 

Figure  4.7:  Amplitude  of  Deflection  Signal  v.v.  Chopper  Frequency 
for  Small  Aspect  Ratio  Reed. 

Deflection  Decrease  With  Increased  Absorber.  The  deflection  signal  of  the  reed  sensor 
went  down  when  the  amount  of  absorber  went  up,  as  shown  below: 
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There  are  several  likely  causes  for  this.  One  possibility  is  that  so  much  of  the  signal  is 
due  to  silicon  absorption  that  the  absorption  of  energy  by  the  sample  results  in  less  heat 
transferred  to  the  overall  reed.  With  the  absorber  on,  the  silicon  itself  absorbs  less,  hence 
less  deflection.  Another  cause  is  that  the  application  of  the  sample  to  the  reed  itself.  The 
sample  may  very  well  be  so  thick  that  the  energy  absorbed  by  the  top  of  the  sample 
transfers  down  through  the  sol-gel  layer  and  to  the  rest  of  the  reed.  A  more  likely  cause 
is  that  the  laser  dye  used  did  not  become  “trapped”  in  the  nanopores,  but  instead  formed 
another  layer  of  its  own  atop  the  sol-gel.  This  addition  of  mass  would  be  expected  to 
decrease  the  deflection  of  the  reed. 
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Mirror.  Early  on  in  the  project,  there  was  concern  that  without  some  sort  of  mirror  or 
reflector  over  the  cantilever  reeds  they  would  not  absorb  enough  energy  to  deflect 
noticeably.  This  is  certainly  the  case,  for  the  monochromator  reduces  the  light  several 
orders  of  magnitude  and  the  fiber  optic  line  by  at  least  six.  Current  mirrors  are  not 
reliable;  made  of  aluminum  foil,  they  can  only  be  re-used  once  or  twice.  Polished  brass 
mirrors  were  constructed,  but  tested  no  better  than  the  aluminum  foil. 

The  mirror  had  no  effect  for  reeds  less  than  1.5  centimeters  in  length.  The  mirrors 
were  proven  effective  for  longer  reeds,  and  enhanced  the  deflection  signal  by  about  five 
times.  This  is  probably  due  to  the  size  of  the  light  beam  and  the  height  of  the  mirror 
After  hitting  the  reed  once,  the  light  is  already  exiting  the  other  end  of  the  mirror,  as  seen 
below.  Lowering  the  mirror  should  solve  this  problem,  but  the  IR  beam  is  currently  so 
wide  that  the  mirror  cannot  be  lowered.  To  contain  the  light  beam  further,  a  wall  can  be 
built  into  the  far  end.  The  wall  would  reflect  the  light  back  to  its  origin  and  should  make 
the  mirror  useful  for  small  reeds. 

Base.  The  reed  has  been  held  on  with  a  variety  of  materials,  including  adhesive  and  a 
clamp.  The  adhesive  reacts  with  the  solvent  for  the  laser  dye,  and  so  this  method  was 
quickly  abandoned.  The  clamp  worked  well  except  for  the  tendency  of  the  jaws  to  snap 
down  on  the  reed  with  too  much  force,  snapping  the  reed.  Also,  the  lips  of  the  clamp 
made  it  hard  to  fit  small  reeds  into  the  AFM.  This  was  overcome  by  filing  the  lips  of  the 
clamp,  which  made  it  of  proper  size,  but  also  made  it  harder  to  clamp  the  reed  to  the  base 
without  snapping  the  reed.  One  way  of  anchoring  the  reed  to  the  base  is  to  insert  the  base 
into  the  AFM  and  then  lay  the  reed  on  top.  Massive  (relative  to  the  reed)  material  with 
similar  dimensions  to  the  base  is  placed  on  top  of  this  assembly.  Small  pieces  of  silicon 
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were  found  to  be  heavy  enough  to  hold  the  reed  in  place.  This  method  is  by  far  the 
easiest  on  the  reed.  It  is  also  easier  to  produce  reproducible  results  as  the  reed  moves  far 
less  with  placement  of  this  weight  than  it  does  with  the  clamping  of  the  reed. 

Position  of  the  Light  Beam.  In  a  perfect  reed,  the  heat  would  be  transmitted  so  quickly 
that  it  would  not  matter  where  the  IR  hit.  This  is  not  the  case,  for  the  diffusion  of  heat 
through  the  first  layer  is  extremely  slow.  Below  is  a  table  showing  the  results  of 
sampling  the  reed  at  three  different  positions.  Position  A  is  closest  to  the  AFM  tip, 
Position  B  is  in  the  middle  of  the  reed,  and  Position  C  is  near  the  base  of  the  reed,  where 
it  is  clamped  down.  This  test  was  run  on  two  different  reeds:  The  first,  a  two-layer  reed 
of  silicon  (5  microns)  and  aluminum  (0.5  microns)  with  measurements  2  centimeter  x  1.5 
mm).  The  other  was  a  reed  of  all  four  layers,  sol-gel  (1  micron),  gold  (2  angstrom), 
silicon  (8  microns),  and  aluminum  (0.5  microns)  measuring  2.25  cm  x  3  mm. 


Figure  4.9:  Deflection  Affected  by  Position  of  Fiber 
In  addition  to  a  decrease  in  signal  as  the  beam  moves  further  away  from  the  base, 

notice  the  marked  decrease  in  signal  that  comes  with  adding  extra  layers.  The  signal 

from  reed  two,  the  Si/Al  reed,  is  what  would  contribute  to  background. 
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Reliability.  The  reed  housing  plays  a  large  role  in  the  repeatability  of  the  measurements. 
There  is  a  small  stopper  that  sits  beneath  the  free  end  of  the  cantilever  reed.  It  was  built 
into  the  reed’s  base  as  a  way  to  keep  down  unwanted  vibrations  (read:  noise)  but  now 
serves  two  purposes:  Reliable  results  are  obtained  only  after  positioning  the  AFM  tip 
over  the  area  of  contact  between  the  cantilever  reed  and  its  stopper.  This  is  because  the 
reed  itself  is  not  stiff  enough  to  hold  up  the  AFM  tip,  and  so  will  continue  to  bend  if  there 
is  not  a  stable  stage-like  presence  beneath  it.  By  taking  care  to  position  the  tip  in  the 
middle  of  the  reed  over  the  stopper,  both  problems  can  be  solved. 

It  was  also  found  important  to  position  the  reed  in  precisely  the  same  place  each 
time.  Figure  4.20  shows  three  data  runs  with  the  same  reed,  a  four-layer  reed  two 
centimeters  in  length.  The  y  axis  signal,  measured  by  the  lock-in  amplifier,  is  the  voltage 
component  of  the  deflection  signal  at  the  phase  indicated  on  the  x  axis.  This  phase  is 
relative  to  the  chopper  frequency,  and  was  taken  from  45°  to  -110  °.  The  intent  of  this 
particular  experiment  was  to  characterize  the  phase  lag  associated  with  this  reed,  but 
complete  measurements  were  not  obtained.  The  first  two,  Run  1  and  Run  2  were  both 
positioned  in  the  same  place  on  the  reed  base.  The  third,  Run  3  hung  about  a  millimeter 
over  the  clamp  end  of  the  base,  effectively  reducing  the  length  of  the  reed.  Notice  the 
extreme  difference  in  signal. 
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Repeatability 
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Figure  4.20:  Component  of  Deflection  Signal  at  Phase  Differences  Relative  to  Chopper 
Frequency.  Shows  the  Importance  of  Positioning  the  Reed  Carefully 


Comparison  of  Model  and  Experimental  Results 

The  laboratory  results  appear  to  validate  the  model,  except  when  important 
assumptions,  such  as  aspect  ratio  of  the  beam,  are  violated.  This  points  to  a  need  for 
more  development  of  the  thermal  problem.  The  reed  sensor  needs  to  be  described  as  a 
plate,  rather  than  a  beam. 

The  model  currently  gives  good  qualitative  analysis;  correctly  predicting 
temperature  profile  and  deflection  mode  shapes,  but  is  not  yet  adapted  for  quantitative 
predictions. 
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5.  Conclusions 


Implications  of  Model  Results 

According  to  the  mathematical  model,  it  is  best  to  operate  with  smaller  reeds. 
Longer  reeds  are  too  sensitive  to  low  frequency  background  noise,  and  shorter  reeds 
produce  too  small  a  deflection  signal.  The  incoming  light  should  be  modulated  near  the 
reed’s  mechanical  frequency.  This  resonance  must  be  found  experimentally,  because 
higher  order  resonance  effects  may  be  present,  depending  on  the  width  of  the  reed. 

The  aluminum  layer  thickness  should  be  increased.  Unfortunately,  current 
fabrication  methods  introduce  significant  curling  of  the  reed  if  the  aluminum  layer  is 
made  any  thicker.  With  improved  construction  techniques,  the  ratio  of  aluminum  to 
silicon  thickness  can  be  brought  to  the  optimum  0.325,  rather  than  the  current  0.0625. 

Implications  of  Experimental  Results 

Experimental  data  validates  this  model  for  use  with  long  thin  beams.  This  model 
is  not  as  effective  for  shorter,  wider  beams  are  easier  to  produce  and  tend  to  better  hold 
the  deposited  sample.  For  this  case,  a  model  based  on  plate  mechanics  rather  than  beam 
mechanics  is  necessary. 

Conclusion 

This  research  was  successful  in  applying  mathematical  models  to  the  reed  sensor. 
The  model  is  applicable  to  long,  thin  reeds  but  does  not  accurately  predict  the  behavior  of 
short,  squat  reeds.  In  lab  experiments  show  that  more  detailed  modeling  is  necessary  to 
describe  these  cases  as  higher  order  effects  become  too  significant  to  continue  modeling 
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the  reed  sensor  as  a  beam.  Modeling  the  reed  sensor  as  a  plate  would  allow  for  the 
development  of  models  that  would  be  valid  over  a  broader  range  of  reed  parameters. 

Suggestions  for  Further  Research 

The  silicon  background  signal  swamps  the  deflection  due  to  any  sample  present 
on  the  reed.  As  mentioned  earlier,  a  decrease  in  deflection  with  increased  absorber  shows 
that  the  reed  sensor  itself  is  more  sensitive  to  light  absorbed  in  the  silicon  layer  than  it  is 
to  light  absorbed  in  the  sample.  This  problem  must  be  eliminated,  perhaps  by  use  of  a 
silicon  filter.  Also,  it  may  be  beneficial  to  consider  alternative  materials  for  this  second 
layer.  For  example,  diamond  is  relatively  transparent  to  IR  light  and  will  not  produce 
such  a  large  background  signal. 

Transmission  of  light  from  the  IR  source  to  the  reed  is  inefficient.  The  optical 
fiber  introduces  a  6  order  of  magnitude  decrease  in  beam  strength.  A  similar  decrease 
results  from  use  of  a  monochromator.  Often,  there  is  not  enough  power  deposited  to 
cause  deflection.  Ideas  to  increase  the  transmitted  fraction  of  light  include  shorter  fibers, 
higher  quality  fiber  tools,  and  light  pipes.  The  first  two  would  bring  obvious  increases  in 
transmitted  fraction.  The  light  pipe,  made  of  highly  reflective  material,  would  deposit  a 
greater  fraction  of  beam  energy  to  the  reed,  but  may  also  spread  the  beam  enough  to 
balance  any  improvement.  Also,  the  relatively  large  diameter  of  such  a  pipe  may  spread 
the  IR  beam  so  large  that  too  little  light  hits  the  cross-sectional  area  of  the  reed. 
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Appendix  A:  Useful  Data 


The  following  table  presents  the  values  that  were  used  in  modeling  the  reed.  The 
numbers  in  parentheses  represent  the  possible  range  for  that  quantity. 

Table  A.  1:  Table  of  Important  Material  Parameters  [1,2] 


Quantity  (symbol) 

Units 

Sol-Gel 

MgAl204) 

Silicon  (Si) 

Aluminum  (Al) 

Thickness  (t) 

(actual  range) 

jLim 

1  d-2) 

8  (4-15) 

.5  (*) 

Density  (p) 

mmm 

3580 

2328 

Elastic  Modulus  (E) 

2.77 

1 

.8 

Thermal  Conductivity  (X) 

W/[m  K] 

14.6 

150 

237 

Thermal  Expansion  (a) 

6.5 

2.6 

23.6 

Heat  Capacity  (Cp) 

J/[kgK] 

837.360 

700 

Thermal  Diffusivity  (D) 

10'b  m2/s 

4.87 

92.05 

96.6 

The  following  table  shows  error  introduced  by  approximating  the  derivative  functions. 
The  error  grew  quickly  for  very  small  reeds,  but  these  reeds  are  outside  of  our  parameter 
space. 


Table  A.2:  Error  in  Derivatives  with  Approximation 


Data  Point 

Baker’s  Values 

Approximated  Values 

%Error 

Root[.03,14,12] 

0.00370472 

0.00370485 

.004 

Root[.0225,8,25] 

0.00279942 

0.00280037 

.03 

Root[.01, 1,500] 

0.000811243 

0.000823453 

1.5 

Root[.004,29,250] 

0.00508951 

0.00402847 

26 
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Appendix  B 

DERIVATION  OF  THE  SOLUTION  TO  THE  THERMAL  PROBLEM 


Fundamental  Equations  and  Conditions. 

Dr  William  Baker,  of  the  Air  Force  Institute  of  Technology,  has  developed 
equations  describing  the  heat  transfer  within  the  reed.  They  are  based  on  the  fundamental 
heat  transfer  equation,  written  for  all  three  layers  of  the  reed.  The  equation  follows  a 
sketch  of  the  reed,  as  seen  from  the  side. 


The  heat  transfer  equation,  written  for  all  three  layers  must  be  solved  to  determine 
the  temperature  profile.  Here  the  concern  is  the  vertical  and  horizontal  directions.  The 
assumption  is  made  that  the  temperature  profile  does  not  vary  with  width.  Further,  it  is 
assumed  that  the  temperature  at  the  left-most  boundary  of  the  reed  is  held  constant  by  the 
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clamp  mechanism  holding  it  in  place.  This  temperature  is  also  the  initial  temperature  of 
the  reed.  The  other  three  boundaries  of  the  reed  are  assumed  insulated.  Convection  and 
radiation  remove  a  negligible  amount  of  heat  from  the  reed. 

The  solutions  for  the  three  layers  are  coupled  together  by  the  additional  conditions 
of  continuity  of  flux  and  temperature  across  the  interfaces.  These  initial  and  boundary 
conditions  are  summarized  as  follows:  For  the  three  insulated  edges,  that  is,  top,  bottom, 
and  right  side: 


The  top  edge  is  insulated: 


du g> 

dy 


(x,h1,t)  =  0 


The  bottom  edge  is  insulated: 


a«(3) 

dy 


( x,-(h2+h3),t)  =  0 


The  right  edge  is  insulated  (requiring  an  equation  for  each  layer): 

a«(1)/r  a«(2)„  x  a«<3\,  x  „ 

— — (L,  y,t)  =  — — ( L,y,t)  =  — — (L,y,t)  =  0 
ox  dx  dx 


The  left  side,  held  at  constant  ambient  temperature,  is  represented  by: 

u (1)  (0,  y,t)  —  u (2)  (0,  y,t)  =  u (3)  (0,  y,  t)  =  U0 


Continuity  of  Temperature  across  the  interface  requires: 

u (1)  (x,0,  t)  =  u (2)  (x,0,  t)  and  u (2)  (x-h2  ,t)  =  u (3)  (x,-h2 ,  t) 


Continuity  of  flux  across  the  interface  further  requires: 

D\~—{xfi,t)  =  D2^ — (*A0  and  D2 ^ — (x,-h2,t)  =  D 3^^-(x,-h  2,t) 

oy  dy  dy  dy 


Finally,  the  reed  starts,  in  a  straight  position,  at  some  initial  temperature: 

u(l)(x,  y,0)  =  UQ 
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Scaling  the  Problem. 

For  ease  of  solution,  the  equation  is  scaled:  lengths  are  scaled  to  the  overall  length 
of  the  reed,  and  time  is  scaled  to  the  period  of  the  modulated  IR  light  beam.  The  equation 
is  re-written  to  reflect  the  real  measure  of  interest:  change  in  temperature.  Assuming  that 
the  reed  starts  at  some  constant  (ambient)  temperature,  the  solution  now  must  be  found 
for  the  deviation  of  that  temperature  with  time. 

If  we  set  x  =  L%,  y  =  Lrj,  t  =  T%  /i,  =  a„  and  define 

u(l)(x,y,t)  =  U0  +  w(%,ri,  r) 


then  the  equations  and  boundary  conditions  become: 

aw(t)(^,77,r) 


dt 

F(T)  =  kTE0(rr) 
dl  =  T^ 


=  rf,V  V°  (|,r7,T)  +  dnF( t) 
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Insulated  edges: 

e  dw^  dw^  dw^ 

-(!,«!, t)  =  ^^-(1,i?,t)  =  — ^  (1,7], t)  =  &-(a2  +«3Xt)  =  0 


dri 


dr] 


Edge  at  constant  ambient  temperature: 

w(1)  ( 0,T],T )  =  W(2)  (0,7],T)  =  W(3)  ( 0,7J,T )  =  0 


Continuity  of  temperature  across  the  interface  gives: 

w(1)  (£  ,0,  t)  =  w(2)  (E,  ,0,  t)  and  w (2)  ,-a2 ,  t)  =  w(3)  (E,  ,-a2 ,  t) 
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Continuity  of  flux  across  the  interface  further  requires: 

and  d2^—^,-a2,r)  =  d3^—(^,-a2,T) 


Finally,  the  reed  starts,  in  a  straight  position,  at  some  initial  temperature: 

°)  =  0 


Assuming  a  Solution. 

To  satisfy  the  boundary  conditions, 

dw(i) 


(1,7J,T)  =  0  and  w(l)(0,77,T)  =  0  for  i  =  1,2,3. 


the  solution  must  have  the  form: 


=  X?(°  (ri,r)Sin(\kZ) 

k=  0 


where 


K 


4=(2k  +  l)-  for  it  =  0,1,2...  . 

2 

This  assumed  solution  is  then  inserted  into  the  initial  equation,  resulting  in 


3w(l) 


0T 

which  becomes. 


k=0 


d<lk 


(«) 


a 


(B.2) 


5t  +  dih2<lkl)  ~di  4*2  f =  8nF(^ 


to  dil2  h 

because  of  orthogonality  of  the  sine  function  over  the  interval  [0,1]. 
Similarly  treated,  the  initial  and  boundary  conditions  become: 
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dqP  dqW 

=  -^-(-(a2  +  a3  ),T)  =  0 

?*)(0,t)  =  ^2)(0,t) 

^2)(-a2,T)  =  ^3)(-a2>T) 


,  &P  a# 

rfi  (0,r)  =  d2  (0,t) 


077 


077 


Mi 

077 


dq™ 

(-a2,T)  =  d3— ^-(hz2,t) 
dr\ 


Laplace  Transform 

The  series  of  qk's  are  Laplace  transformed  with  respect  to  t.  Defining, 

G*°  OM)  =  L^°  (t?,T)j=  Jo°°e_  JT^°  (J/,T)</T , 
the  differential  equations  become 


a2ft(0  2gM7(T) 

'  sij2  4 


Or 


where 


C+rfAba® 

dr/ 


-2SaF(T) 


F(T)  =  j~e-"F(T)dT. 


A  solution  is  easily  obtained  in  the  form: 

W 

where 

Hi=y[s  +  d~A l  . 
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The  unknown  coefficients,  A,  and  Bj,  are  determined  by  applying  the  transformed 


boundary  conditions, 

.  30®  ae<3) 

i.  -^(a1,s)  =  -^—(-(a2+a3),s)  =  0 
»•  Gi1)(O,5)  =  0f>(O,5) 

30®  ao[2) 

in.  dx  — MO,  j)  =  d2  -¥*-(0  ,s) 
orj  dr\ 

•  ,  30*®  aoJ3) 

*v-  =  d3  -^-(-a2,s) 

v.  Qf\-a2,s)  =  Qf\-a2,s) 

Defining, 


7i  = 


■id, 


and  employing  boundary  condition  i  results  in: 

0?»s)  =  A  (s)Cosh(-jL  ( r ] -ax))  +  — 

y]dl  Xk!i{ 

—n  —  r) 

Qjc2)(T],s)  =  A2(s)e ^  +  B2 (s)e  ^ 

Q(k]  07 » ■*)  =  A  (s)Cosh(-j=  (j]+a2+  a3 )) 

V^3 


Applying  boundary  conditions  ii  and  iii,  along  with  some  algebra,  yield: 

M» = - ^4  Co^ftr, ) -^-WsinKw, ) 

2  2  ^  ^2VJ2 

2  2  l,  ^2V«2 

after  defining 


#*(*)  = 


2F(£) 
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Doing  the  same  for  conditions  iv  and  v  results  in: 


£ 


Cosh(n3y3 )  +  —J^~Sinh(ju3y3 ) 
M2  V« 2 

/ 

Cosh(n3y3 )  -  —  \^j-Sinh(ji3y3 ) 

M2  \d2 


If 


ci  =  Coshlpj-)  and  S;  =  --^-Sinhiju^) , 

M2\d2 


then  the  previous  results  (combined  and  rearranged)  become 

M*)  =  \Hk  (s)  +  ^  Ai  )(Q  -Sy)  =  As (5)(C3  +  ^  ) 

B2(s) =  —H  k  (s)+-Ms)(Ci+S1)  =  ^e~M2Y2A3(s)(C3-S3) 


And  can  be  manipulated  to  obtain: 


Q  ~(C2C3+S2S3)  Aj(>y) 

Si  (S2C3+C2S3)  _  _A3(.y)_ 

0 

The  coefficients  Ai  and  A3  are  easily  solved  as, 

AC5)  =  ”'(‘^2^3  ^2*^3 ) 

A3(s)  =  j±S, 

where 

D  =  Cl  ( S2C3  +  C2S3)  +  Sj  (C2C3  +  S2S3 ) 
and  substituted  to  solve  for  A2  and  B2: 
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AI(S)=Ie',^^£L(C3+S3) 

B2(s)^e~^2lh^lic]-S,) 


Note  the  distinction  between  s,  transformed  time,  and  S,-,  a  function  written  for 
convenience.  Finally,  the  solution  for  each  layer,  in  the  transform  domain  is: 

Qk>  (V,s)  =  -~(S2C3  +  C2S3  )Cosh(m  (-?=  -7i))+ Hk 
u  yld  l 


Qf]  (*M)  =  Cosh(-j==  (7 n+a2+  a3 )) 

D  Vd3 


i \ = 


D 


1  M2(72+-n=)  -^2(Y2+~T=) 

±C3(e  ^2  +e 


)  - 


1  »2(ri+-j=)  -n2(Y2+-rr) 

+  -S3(e  ^2  -e 

2 


D 


C3Coshtx2  (y2  +  -j=)  +  S3Sinhfl2  (y2  +  -^L) 
V“2  v^2 


Careful  examination  shows  there  are  no  branch  points,  and  further,  that  if  z  is  a 
root,  then  z  is  as  well. 


Inverting  the  Laplace  Transform. 

Directly  inverting  the  Laplace  transform,  via  calculus  of  residues,  requires  finding 
the  roots  of 

D(s)  =  fi2(s)D{s)  =  n2{s)(ClS2C3  +  ClC2S3+SxC2C3+SlS2S3). 
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For  the  parameters  of  this  problem,  D(s)  turns  out  to  have  a  single  root  for  each 
k;  lying  along  the  negative  real  axis. 


Defining, 


f^i  —  J 


(O 


(0 


1(0 


(0 


zk  ~  ~°k 


if  di^l~ak>0 
if  dfXl-ak<  0 


The  inverse  Laplace  transform  of  a  representative  piece  of  the  solutions  becomes 

L~l  Cosh(Hi  (. s)fi  (77))  j 

1  (c+i°°„st  Ms) 


where 


_  1  fc+i 

9 777  Jc-I« 


2mJc-ioo  £>(5.) 


Cosh(fli(s)p(rj))ds 


=  ^^CoshiHi  (5)78(77))  j 


=  e~akt  A(~a k)!1!^ k)  \Cosh((o[l) (5)78(77))  if  J. A*  -ck>0 

D'i-^k)  Cos(co£\s)P(ri))  if  di7?k-ck>  0 


D'(s)^(5(s))  =  -f(fi2(s)D(s)) 
as  ds 


Within  the  parameter  space  of  this  problem,  these  roots  and  functions  evaluate  to 


zk  -  “Ok 
fii(Zk)  =  ia>k) 
/*2 (zk)  =  0)jc2)  . 
fh( zk)=co 
At  =  D'(-ok ) 


This  last  term  reflects  the  derivative  function  being  real  at  each  root. 

So  that  the  convolution  theorem  may  later  be  used,  the  Qk’s  are  written  as 
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G®  (ij,  J )  =  (S2C3  +  c2s3  xw(ft 

( 

Gk2)  <n> *)  =  ~r^  C3Coshii2 (72  +  -J=)  +  S3Sinhfi2 (y2  +  -^L) 

V“2  yjd  2 

th 


KD 


Gf  (?],  =  — f-  Cosh(-j3=  (xi +a2+  a3 )) 

kd  Vrf3 


Now,  after  defining 


gf(77,^)  =  L-1{Gf07,5)J 

n\  1  /Y)^ 

gl)(V>s)  =  -e~°kt  -  Cos(n= (7? + 1 ®i »  * 

|  (0(f)Sinh((of)y2  )Cosh((of)y3 )  +  J^co(k3)Cosh(o)(k2)y2  )Sinh(cof)y3 )  j 


\d2  (o(pDk 


Cosh((of^y3 )Cosh(^j=  (rj + y23))  +  Sinh(cof)y3 )Sinh(°^L (77  +  y2  )) 

2  V  rf2  ©f ;  yd2 


,(3) 


(3),  v  —Oirt  [d7  Sinh((o9^y,)  a)[ 

g[  \-n,s)  =  e-°kt  ~ - Cosh(-j=(T]  +a2+a3)), 

\d2  (0\iDk  yds 


the  <yf (77,0  are  determined  by  convoluting  gf  (77,7)  with  — F(t).  To  capture 


the  time  dependence  of  the  gf  (77,  t) ,  Itft)  is  defined  as 


Ik(t)  =  j‘e-^F(t)dr 


Incoming  Energy. 

It  is  time  to  consider  the  absorbed  energy,  F(t).  This  energy,  originating  from  light 
modulated  by  a  chopper,  takes  the  form  of: 


F(T)  =  kE0(x)  =  k\ 


Ir\  0<  5s  tr\ 


0  tn<—  <1 


(B.3) 


Here  tQ  represents  the  fraction  of  a  period,  T,  that  the  beam  is  turned  on.  Only  a 
portion,  k,  of  the  incident  energy,  E0,  is  deposited  in  the  sol  gel  layer. 

Up  to,  and  including,  one  period,  the  solution  for  the  time  dependant  part  of  the 
problem  looks  like 


\l-eakt  0<T  <t0 

\(e°ktQ  -\)e~Gkt  t0<t<T 


Of  course,  the  following  form  is  needed  for  t  over  many  (more  than  1)  periods: 


hit) 


°k 


(eCT*t0  -1) 
{e°k  -1) 


0 e°kP  -\)-e°kP 


kl0  f  1  0  <  T  <t0  (B-4) 

+  CT^  t0<T< 


where  p  is  the  number  of  periods  completed.  The  parameters,  k,  and  Io  can  be  determined 
experimentally  as  described  in  Chapter  3. 


Min  and  Max. 

Casual  analysis  of  the  time  dependant  portion  of  the  solution  indicates  that  the 
greatest  change  in  temperature  will  occur  when  to,  the  portion  of  the  period  that  the  beam 
is  turned  on,  is  equal  to  0.5.  The  chopper  used,  by  design,  is  set  up  such  that  to  =  0.5. 
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Summary  of  Solution. 

The  final  solution  for  each  layer  is  comprised  of  a  series  of  solutions.  Each  series 
term,  k,  includes  a  dynamic  description  of  the  energy  deposition, 


(e°k  -1) 


(e°kP  -\)-e°kP 


kin 


1 


0<~<tn 


.  < l<i 

0  T 


modulated  in  the  horizontal  (x)  direction  by  (f>k  (£)  =  —Sin{hkt;)  . 

K 

The  energy  attenuates  as  it  passes  through  the  depth  (y)  of  the  reed.  This  behavior 
is  shown,  for  each  layer  i,  by  y/‘k  (rj) : 


i  (o? 

YkOl)-  Cos(—j=  {r\ + ax )) 

V“i 

2  <o[2)  m,(2) 

¥kOl)  =  AkCosh(-^=('n  +  a2))  +  BkSinh(-^(T1  +  a2)) 

V“2  V^2 


(B.5) 


where 


Ak  =  Cosh(a>^y3) 

fV\ 

— -qT  Sinh(G)jpY3 ) 

cop 


0) 


(3) 


Yk  (7)  =  Cosh(—jL=  0 n+a2+  cx3 )) . 
■yjdj 


Finally,  the  solution  to  the  thermal  behavior  for  a  three-layer  reed  can  be  written: 


=  Xc*0/*  W*  £)yT  (n) 


k=  0 


(B.7) 


B.12 


where 


r(D  - 


-2 


(0>P)2Dk 


W2)Sinh((0(2)y2  )Cosh((0ik)Y3)  +  .f^-0)^>Cosh(o)j2)y2  )Sinh((of)y 3 )] 


and 


cf=cf  = 


[d7  Sinh((O^Yi) 
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Appendix  C:  Mathematical  Model 


Thermal  Solution 
Approximation  of  the  Code 
Application  to  the  Reed  Sensor 
Phase  Exploration 

Resonance  and  Optimum  Thickness  Exploration 

■  Reed  Parameters 

First,  the  important  thermal  parameters  are  defined:  thickness  of  each  layer,  length, 
frequency  and  period  of  the  chopper,  thermal  diffusivity 
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(*  height  of  layers  *) 
hi  =  1*10A  (-6); 
h2  =  8  *  10  A  (-6)  ; 
h3  =  . 5 *  10 A  (-6); 

(★length  of  reed, 

assuming  Aspect  (L/w)  ratio  of  at  least  10  *) 
L  =  .0225; 

(★Modulation  frequency#  in  Hertz*) 
ChopperSpeed  =8; 

(★Timescale,  in  seconds  ★) 

T  =  N[1  / ChopperSpeed] ; 

(*  alphas,  scaled  height,  unitless  *) 
alphal  =  hi  /  L; 
alpha2  =  h2  /  L; 
alpha3  =  h3  /  L; 


(*  thermal  diffusivities,  in  mA2/s  *) 
D1  =  4 . 87  *  10 A  (-6); 

D2  =  92.05*  10  A  (-6)  ; 

D3  =  96.6*  10A  (-6)  ; 


(*  scaled  diffusion,  unitless  *) 
ddl  =  (T*D1)  /LA2; 
dd2  =  (T*D2)  /LA2; 
dd3  =  (T*D3)  /LA2; 


(★  "diffusion  speed"  scale. 


garni  = 


gam2  = 


gam3  = 


alphal 
-  • 

Vddl 

alpha2 

Vdd2 

alpha3 

Vdd3 


unitless  ★) 


■  The  Root  Solving  Block 

Here  the  relationships  are  defined  and  the  root  terms  are  calculated.  The  derivative 
terms  are  also  calculated. 

A*  =  (2k+l)*f  and  "lams"  =  A*2 
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Px  ^ 

lams[k_]  :=  ( (2*k+ 1)  * - J  ; 


Mi  =  yfs  +  d^f 


itml[s_,  k_] 
mu2  [s_,  k_] 
mu3  [s_,  k_] 


Cj  =  Cosh  {Hi  ji) 


:  =  V  s  +  ddl  *  lams  [k]  ; 
:=  V s  +  dd2  *  lams [k]  ; 
:=  V  s  +  dd3  *  lams  [k]  ; 


ccl[s_,  k_] 
cc2  [s_,  k_] 
cc3[s_,  k_] 

si~  ^-Sinh (jiiji) 


:=  Cosh[inul[s,  k]  *gaml]; 
:=  Cosh[mu2  [s,  k]  *  gam2]  ; 
:=  Cosh[mu3  [s,  k]  *gam3]; 


ssl  [  s _ ,  k_]  !  = 


V 


ddl 

dd2 


mul[s,  k] 

*  — — - —  *  Sinh[mul[S/  k]  *  garni]  ; 

mu2 [Qf  k] 


ss2[s_,  k_]  :=  Sinh[imi2 [s,  k]  *gam2];  taken  out*) 

ss3  [  s_,  k_]  :  = 


mu3[s,  k] 

* - 

mu2[s,  k] 


*  Sinh[mu3  [s,  k]  *gam3]; 


The  function  we  want  roots  to:  D  =  /x2( s)D(s)=  cx  c3  s2  m  +  cx  c2  53  fi2  +  c2  c3  s{  /i2 
+  Si  s3  s2  n2 


M[s_,  k_]  ^cclts,  k]  *cc3[s,  k]  *ss2[s,  k]  *mu2[s,  k] 

+  ccl[s,  k]  *cc2[s,  k]  *ss3[s,  k]  *mu2[s,  k] 

+  cc2[s,  k]  *  cc3  [s,  k]  *ssl[s,  k]  *mu2[s,  k] 

+  ssl [s,  k]  *  ss2  [s,  k]  *  ss3  [s,  k]  *mu2[s#  k] 

Here  the  root  findng  block  will  find  the  roots,  assign  them  to  the  values  <x*,  and  also 
assign  the  corresponding  values  to 

<3  ■§■  s  ( zk )  /  Mi  ( zjc ) ,  M2  (Zk) ,  and  M3  ( zk )  •  Mathematica  uses  Newton's 
method  to  search  for  the  root. 


Set  k  max  before  running  it.  Also,  note  that  the  table  which  lists  all  of  these  values  is 
not  printed  out,  you  can  easily  have  it  print  by  removing  the  comment  markers  "(*" 
and  "*)"  which  bracket  it. 
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info  =  {{  "k"  ,  “A*",  »-d2Ak2",  "zk", 

"4^(zO%  n^i(zk)",  n^(zk)",  "/i3  (zk) "}}; 
as  J  J 

kxnax  =  100; 

Do[ak  =  Vlams[k]  ;  bk  =  -dd2*  lams[k]  ; 
ck=  z  /•  FindRoot  [Dd[z,  k]  ==  0,  {z,  .99*bk}]  ; 
slg[k]  =  Abs[ck]; 
dk  =  (D[Dd[s,  k]  #  s]  / .  s  ->  ck)  ; 
derv[k]  =  dk; 
ek  =  ittul  [ck,  k]  ; 
wl[k]  =  Im[ek]  ; 
fk  =  mu2  [ck,  k]  ; 
w2  [k]  s  Re[fk]  ; 
gk  =  imi3  [ck,  k] ; 
w3[k]  =  Re[gk] ; 

AppendTo[ 

infO|  {k,  ak,  bk,  ck,  dk,  6k,  £k,  gk} } ;  ,  {k,  0,  knax}  J 
(*NuniberForm[TableForm[info,TableSpacing->l]  ,7]  *) 

Now  the  series  solution: 


The  Constant  Functions  "ck" 


ckl[k_]  := 
ckl[k]  = 


(wl[k]  )2  *derv[k] 


w2[k]  *Sinh[w2[k]  ★  gam2]  * 


Cosh[w3[k]  *gam3]  + 


dd3 

dd2 


w3  [k]  * 


Cosh[w2[k]  *gam2]  *Sinh[w3[k]  *gam3] 


c*2[k_]  .  Sin[wl[k]  .garni] 


dd2  wl[k]  *derv[k] 


=»[*_,  ,.=«m  =  2.,/  221 ,  .gi] 


dd2  wl[k]  *derv[k] 


■  The  Time-Energy  Behavior  "I" 

This  function  describes  the  deposition  of  energy  into  the  layers.  The  expected 
difference  between  the  layers  comes  from  the  modlation  of  I  by  the  other  terms. 
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kk  =  1; 

(♦this  k  represents  an  absorption  of  the  top  layer*) 
Io  =  1  *  10 A  - 6;  (*  this  represents 

the  initial  intensity  of  the  incoming  beam. 

The  fibre  itself  effects  a  10~6  reduction  *) 
to  =  . 5  (*  fraction  of  period  that  signal  is  on, 
this  cannot  change  and  is  due  to  chopper  setup  *) ; 
P[t_]  :  =  p[t]  =  IntegerPart [t] 


Ii [t_,  k_]  :=  Ii[t,  k]  =  kk*I°  *E-al°W*t* 

sig[k] 

(fEsig^*to  -  \ 

- —  *  (E8iff[k]*p[t]  -  1}  - Eaig[k]*p[t]  + 

(E»isrW  _  1)  '  ) 


kk*  Io 

—  —  *If[t-p[t]  <=to,  1, 
sxg[k] 


E-siB[k]*((t-p[t])-to)]; 


■  The  X  Behavior,  <p,  or  "Thi" 

This  series  behavior  was  assumed  into  the  solution,  and  accounts  for  the  thermal  sink 
on  one  end  and  an  insulated  condition  on  the  other  end. 

0k  (<T)  =  y-  Sine  (Ak  £) 

Ale 

lam[k_]  :=lam[k]  =  ((2*k+l)  *-^j; 

thi[?_,  k_]  :  =  thi [£,  k]  =  - *Sin[lam[k]  *fl; 

lam[k] 


■  The  Y  Behavior,  iff,  or  "Psi" 

Relation  for  the  thermal  behavior  in  the  vertical  direction. 
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Aa[k_]  :=  Aa[k]  =  Cosh[w3  [k]  *  gam3]  ; 


Bb[k_]  :=  Bb[k] 


dd3  w3  [k] 

* - *  Sinh  [w3  [k]  *  gam3]  ; 


dd2  w2[k] 
Psil[?7_,  k_]  :=  Psil[?7,  k]  =  Cos[wl[k]  * 
Psi2  [r)_,  k_]  s  = 

Psi2[r7,  k]  =  Aa[k]  *Cosh[w2[k]  * 


l  Vddl 


gaml 


Bb[k]  *Sinh[w2[k]  * 
Psi3[77_#  k_]  :  = 


{  Vdd2 


V  dd2 

Sam2]  ]  ; 


+  gaxn2 


,  w3  [k] 


Psi3 [rj,  k]  =  Cosh[ - *  (r]  +  alpha2  +  alpha3)  1  ; 

Vdd3 


]’ 


■  The  Final  Equation:  Putting  it  all  together 

WW1[£_,  r\_,  z_]  :=Wwl[£,  rj,  r]  = 

kmax 

2  ckl [k]  *  Xi[r,  k]  *thi[£,  k]  *Psil[r7,  k]  ; 

k=0 

Ww2 ry_#  r_]  s=WW2[S#  rj,  z]  = 

kmax 

ck2[k]  *Xi[z:#  k]  *  thi[^,  k]  *Psi2[r?,  k]  ; 

k=0 

Ww3[S_#  z_]  ;=WW3[ g,  r\,  z]  = 

kmax 

^  ck3  [k]  *  Xi[r,  k]  *  thi[ff#  k]  *Psi3[?7/  k]  ; 

k=0 


■  Sample  Calculations: 

Here,  I  present  some  sample  calculations  using  series  terms  up  to  and  including  100. 
The  data  points  seem  very  small,  this  is  because  of  the  reduction  of  light  intensity  as 
it  travels  through  the  beam.  Dr  Li  calculated  this  to  be  a  10"6drop: 

These  next  few  blocks  are  groups  of  six  data  points.  They  take  place  at  various  points 
on  the  x  axis,  as  marked-  at  one  half  of  the  first  period  (so  the  beam  is  just  shutting 
off)  and  are  located  at  the  boundaries. 

The  first  point  is  on  the  top  of  the  sol  gel  layer. 

The  next  two  printouts  are  the  boundary  between  the  silicon  and  sol  gel  layers  (both 
outputs  should  match,  but  they  don't).  The  second  set  of  points  represents  teh 
boundary  between  the  silicon  and  aluminum  layers  and  again  should  match. 
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Finally,  the  last  data  point  printed  in  each  set  is  a  point  on  the  very  bottom  of  the 
aluminum  layer. 

At  the  end  of  the  reed,  that  is  x=L 

WWl[l,  alpha 1,  2000.5] 

WW1[1,  0,  2000.5] 

Wto2[l,  0,  2000.5] 

WW2[1,  -alpha2,  2000.5] 

WW3[1,  -alpha2,  2000.5] 

Ww3  [1,  -  (alpha2  +  alpha3) ,  2000 . 5] 

1 . 2955  x  10”6  +  0.1 

1 . 2955  x  10~6  +  0.1 

1 . 2955  x  10~6  +  0.1 

1 . 2955  x  10~6  +  0.1 

1 .2955  x  10~6  +  0.1 

1 . 2955  x  10“6  +  0.1 

Wtnrl[l,  alphal,  2000] 

WW1[1,  0,  2000] 

Wto2  [1,  0,  2000] 

Wm2[1#  -alpha2,  2000] 

Vtw3  [1,  -alpha2,  2000] 

Wto3  [1#  -  (alpha2  +  alpha3)  ,  2000] 

1 . 26919  x  10“6  +  0.1 

1 . 26919  x  10"6  +  0.1 

1 . 26919  x  10“6  +  0.1 

1 .26919  x  10"6  +  0.1 

1 . 26919  x  10"6  +  0.1 

1 .26919  x  10“6  +  0.1 
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x  =  .OIL,  that  is,  near  the  beginning  of  the  reed 

Wtol[.01,  alphal,  2000.5] 

Wtfflt.Ol,  0#  2000.5] 

WW2[.01,  0,  2000.5] 

Wto2 [ .01,  -alpha2,  2000.5] 

WW3 [ . 01,  -alpha2,  2000.5] 

Ww3 [ .01,  -  (alpha2  +  alpha3) ,  2000.5] 

2 .76281  x  10~8  +  0.1 
2 .7628  x  10"®  +  0.1 
2 .7628  x  10~8  +  0.1 

2 .76279  x  10"®  +  0.1 

2 .76279  x  10~®  +  0.1 

2 . 76279  x  10"®  +  0.1 

x  =  .OIL,  again,  this  time  at  the  beginning  of  a  period,  instead  of  half  way  through. 
Notice  that  the  layers  are  all  the  same. 

Wtol[.01,  alphal,  2000] 

Wwl[ .01,  0,  2000] 

WW2 [ . 01,  0,  2000] 

Ww2[ .01,  -alpha2,  2000] 

WW3 [ .01,  -alpha2,  2000] 

Wto3  [  .01,  -  (alpha2  +  alpha3)  ,  2000] 

2 . 34155  x  10'®  +  0.1 

2 . 34155  x  10'®  +  0.1 

2 .34155  x  10'®  +  0.1 

2 . 34155  x  10”8  +  0.1 
2.34155x10"®  +  0.1 
2.34155x10'®  +  0.1 
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The  following  graph  is  temperature  vs  scaled  time,  for  the  top  layer: 

Plot  [Abs [Wwl [1,  0,  t]  ]  ,  {t,  0,  30}]; 


The  Below  Graph  is  Temperature  Fluctuation  vs  Depth,  combined  for  all  three  layers, 
The  first  plot  is  at  the  beginning  of  the  2001th  period,  steady  state.  The  second  is 
halfway  through  the  2001th  period.  The  gridlines  represent  the  edges  of  each  layer. 
Note  that  they  both  have  the  same  shape.  At  the  beginning  of  the  period  has 
magnitude  about  1 .269  *  10-6  Kelvin,  and  half  a  period  later,  the  temperature  in  the 
layers  has  jumped  up  to  the  1.295  *  10~6  Kelvin  range. 
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a  =  Plot [Abs [WWl [  1,  y,  2000.5]],  {y,  0,  alphal}]; 
b  =  Plot [Abs [ VJW2 [  1,  y,  2000.5]],  {y,  -alpha2,  0}]; 
c  =  Plot  [Abs  [WW3[1,  y,  2000.5]], 

{y,  -  (alpha2  +  alpha3)  ,  -alpha2}]; 
w  =  Show  [a,  b,  c. 

Frame ->  True,  FrameLabel  ->  {"Depth"  ,  "  AT  "}, 

FrameTicks ->  {{{alphal,  "ti"},  0,  {  _alpha2  ,  »_.^lii}# 

2  2 

{-  (alpha2  +  alpha3) ,  "- (t2+t3)  "}}, 

{.0276281,  .0276279}},  GridLines -> 

{  {alphal,  0,  -alpha2,  -  (alpha2  +  alpha3)  } ,  None}  , 
PlotLabel  ->  "AT  ,  Steady  State,  t=0  to"] ; 


AT  ,  Steady  State,  t=0  to 
t2 

-(t2+t3)  2  o  fcl 


2 

Depth 
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a  =  Plot [Abs [Wtol [  1,  y,  2000]],  {y,  0,  alphal}]; 
b  =  Plot  [Abs [VJW2  [  1,  y,  2000]],  {y,  -alpha2,  0}]; 
c  =  Plot  [Abs  [Ww3[l,  y,  2000]], 

{y,  -  (alpha2  +  alpha3) ,  -alpha2}  ]  ; 
w  =  Show  [a,  b,  c. 

Frame  ->  True,  FrameLabel  ->  { "Depth"  ,  ■  AT  "}, 

-  alpha2  to 

FrameTicks ->  {{{alphal,  "ti"},  0,  { - ,  " - "}, 

2  2 

{- (alpha2  +  alpha3) ,  (t2+t3)  "}}, 

{.0276281,  .0276279}},  GridLines -> 

{  {alphal,  0,  -alpha2,  -  (alpha2  +  alpha3)  }  ,  None} , 
PlotLabel ->  "AT  ,  Steady  State,  t=  .5  to"]; 


AT  ,  Steady  State,  t=  .5  to 


-(t2  +  t3)  2  0  fcl 


2 

Depth 


Approximation 

Here  I  make  some  approximations  to  Dr  Baker's  thermal  code.  I  ran  the  thermal  code 
for  varying  lengths  and  frequencies  and  fit  the  data  as  shown  below: 

■  Frequency 

■  Roots  as  a  function  of  w,  for  k=0 

Here  I  create  a  List  of  the  Frequencies  and  their  Root  Values,  obtained  in  another  file 
using  Dr  Baker's  Thermal  Code  These  are  all  for  the  first  root,  that  is,  k=0: 


C.ll 


rootO =  {{1,  -.4050890),  {5,  -.08101618), 

{8,  -0.05063511),  {15,  -0.02700539), 

{30,  -0.0135027),  {55,  -0.0073651078), 

{80,  -0.005063511130),  {100,  -0.004050809), 
{130,  -0.003116007),  {155,  -0.002613425}); 

a  =  ListPlot [rootO,  PlotStyle  ->  PointSize [ . 02]  ] ; 


Note  that  the  magnitude  of  the  roots  approach  zero  as  the  frequencies  approach 
infinity 


1 

Fit  [rootO,  {  — ,  1,  w},  w] 

9.34387  xlCT7  -  -•405°9  -  7 . 48072 x 10'9  w 
w 

Testing  showed  only  the  middle  term  to  be  signficant  (small  coefficients  are 
significant  in  other  fits) 
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a  =  ListPlot [ root 0 ,  PlotStyle  ->  PointSize[ . 02]  ]  ; 
b=  Plot [.405089  /-x,  {x,  0,  155}]; 

Show[a,  b] ; 


■  Derivative  as  a  function  of  w  for  k  =  0  (pretty  much  constant  all  the  way  through 
k  =  25) 

Here  I  create  a  List  of  the  Frequencies  and  their  Derivative  Values,  11  for  k=0: 

dervO  =  { {1,  0.0009901744},  (5,  0.002214097}, 

{8,  0.002800636},  {15,  0.003834929}, 

{30,  0.005423409},  {55,  0.00734333}, 

{80,  0.008856389}  ,  {100,  0.009901744}, 

{130,  0.01128973},  {155,  0.01232757}}; 

Fit  [dervO,  {l,  Vw,  w},  w] 

-4.76512  x  1CT10  +  0 . 000990175  Vw  -  6.09522  x  10~12  w 


a  =  ListPlot [dervO,  PlotStyle  ->  PointSize[ .02] ] ; 
b  =  Plot [0.0009901744  Vx#  {x,  0,  155}]; 

Show[a,  b] ; 
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■  Length  (all  for  k=0) 

■  Roots 


rootO  =  {{.01,  -2.050722},  {.005,  -8.202887}, 
{.0075,  -3.645728},  {.0125,  -1.312462}, 
{.015,  -0.911432},  {.0175,  -0.6696235}, 
{.02,  -0.5126805},  {.0225,  -0.4050809}, 
{.025,  -0.3281155}}; 

Fit  [rootO,  {-££-}#  L] 


0.000205072 


a  =  ListPlot [rootO, 

PlotStyle  ->  Point  Si  ze[  .02]  ,  AxesOrigin  -  >  {0,  0}]  ; 
,  r  -0001 

b  =  Plot f - *  (-2.050722), 

L  (x2) 

{x,  .004,  .025},  AxesOrigin ->  {0,  0}]; 

Show[a,  b,  AxesOrigin ->  {0,  0}]; 


-4 


-6 


0.005  0.01  0.015  #  0^)2  •  0.«25 


Here  the  fit  is  Z(L)  =  .0001  where  Zo  =  Z(L=.01),  calcuated  by  the  root  solver 
used  in  the  thrmal  code. 


■  Derivatives-  doesn't  vary  much  with  L,  there  is  tiny  (7th  decimal)  variance,  but 
for  the  most  part  all  the  deriviatives  equal  the  derivative  at  L=.01 

dervO  =  {{.01,  0.0009901742}, 

{.005,  0.0009901737},  {.0075,  0.0009901741}, 

{.0125,  0.0009901743},  {.015,  0.0009901744}, 

{.0175,  0.0009901744},  {.02,  0.0009901744}, 

{.0225,  .00099017440809},  {.025,  0.0009901744}}; 
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Fit[dervO,  {l,  ^pL ,  l},  l] 

0.000990172  +  4.08705x10'®  -  1.47073  xlO'7  L 

Li st Plot [dervO,  PlotStyle  ->  PointSize [ . 01]  ] ; 


0.000990174 

•  •  •  •  • 

0.000990174 

• 

0.000990174 

0.000930174 

0.000990174 

0.000990174 

0.000990174 

0.005 


0.015 


0.02 


0.025 


■  Roots  Vary  With  K 

The  following  roots  were  calculated  in  the  file  "varying.nb"  for  L=. 01  with  w=  1. 

roots  = -1*  {2.050722,  18.45649,  51.26802, 

100.4853,  166.1082,  248.1366,  346.5706, 

461.41,  592.6545,  740.3055,  904.3608,  1084.821, 
1281.686,  1494.955,  1724.629,  1970.708,  2233.19, 
2512.075,  2807.365,  3119.057,  3447.153,  3791.651, 
4152.551,  4529.853,  4923.557,  5333.662}; 
a  =  ListPlot [roots] ; 
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Pit  [roots,  {1,  k,  k2},  k] 

-2.02072  +  8.19207  k-  8.20216  k2 

rootf it [k__]  :=  -2.02072  +  8.19207*  (k)  -8.20216*  (k)2; 
b  =  Plot [ root f it [k],  {k,  0,  26}]; 

Show[a,  b]  ; 


Now  to  see  how  well  the  fit  holds  for  further  values.  The  highest  k  I've  used  is  500. 
The  thermal  code  generated:  k=  500,  zk=  -2.007824  *  106.  This  fit  gives  a  value  of 
-2.04645xl06.  Also,  way  out  here  the  Derivative  function  has  changed  from 
0.0009901742  to  .0008234526,  but  by  this  time,  the  root  is  -70,000  so  the  error 
shouldn't  make  much  of  a  difference! 
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■  Combine  the  three  effects: 


■  Roots  (All  Cases  Okay  except  the  Very  small  (.004): 

I  expect  that  Z(L,w,k)=  where  Zo(k) 

=-2.02072  +  8.19207  *  k  -  8.20216  *  k2 

Here  the  fit  is  Z(L)  =  .0001  where  Zo  =  Z(L=.01),  calcuated  by  the  root  solver  for 
each  reed.  Note  that  Mathematica  calls  Zo[0]  Zo[l]  hence  the  use  of  +1  terms 


Zo [k_]  := -2.02072  + 8.19207*  (k+1)  -8.20216*  (k+l)2j 
.0001*  Zo[k] 


Zz[L_,  w_,  k_] 


U*I? 


Here  is  How  my  Approximation  Performed  for  Calculating  the  Roots: 


Data  Point  Actual  Values 
Zz[.01,l,500]:  -2.007824*  106 

Zz[.01,8,25]:  -666.708 

Zz[.0225,8,25]:  -131.7003 
Zz[.03,14,12]  -10.1722 

Zz[.004, 29,250]:  -106512. 
Zz[.225,8,25]:  -1.31702 


Calculated  Values  %  Error 


-2 . 04645  x  106 

1.9% 

-666.711 

.0005% 

-131.696 

.0032% 

-10.1721 

.001% 

-110041. 

3.3% 

-1.31696 

.004% 

■  Derivatives-  All  Cases  Good  Except  the  Very  Small  (.004) 

It  is  seen  that  the  deriviatives  do  not  change  much  with  "L"  and  are  essentially 
constant  once  calculated  for  a  particular  reed  at  a  particular  frequency.  They  do  vary 
very  slightly  with  k,  a  slight  difference  here  which  causes  catastrophic  differences  in 
the  actual  temperature  calculation! 

Dd[L_,  w_,  k_]  :  = 

(0 . 0009902715  -  3 . 328841  *  10‘9  *  IntegerPart  [  k]  - 
7.094576316*  10"10  *  (IntegerPart [k] ) 2)  * 

Here  is  How  my  Approximation  Performed  for  Calculating  the  Roots. 


Data  Point 

Fitted  Model 

Actual  Values 

%  Error 

Dd[.01, 1,500]: 

0.000811243 

0.000823453 

1.5% 

Dd[.0 1,8,25]: 

0.00279942 

0.00279928 

.005% 

Dd[.0225,8,25]: 

0.00279942 

0.00280037 

.03% 

Dd[.03,14,12] 

0.00370472 

0.00370485 

.004% 

Dd[.004, 29,250]: 

0.00508951 

0.00402847 

26% 

Dd[.225,8,25]: 

0.00279942 

0.00280063 

.04% 
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Here,  the  values  in  the  list  "dervs"  have  been  calculated  in  the  file  "varying.nb"  for 
L.01  and  T  =  1/  1  Hz  and  were  imported  in  without  being  displayed.  The  fit  is: 

derv  =  Flatten  [dervs]  ; 

Derv[k_]  :=  Fit  [derv,  {1,  x,  x2},  x]  /.x->k 
Derv[k] 

0.000990272  -  3 . 32884  x  1CT9  k  -  7 . 09458  x  1CT10  k2 

The  fit  was  using  a  data  set  up  to  250.  Now  let's  check  predictability: 

Derv[k_]  s  = 

0.0009902715-  3.328841*  10'9  *  IntegerPart  [  k]  - 
7.094576316*  10“10 *  (IntegerPart [k] ) 2; 
b  =  Plot  [Derv [k] ,  {k,  0,  500}]; 

Show[a,  b] ; 


Application 

Now  to  Apply  Dr  Baker's  Solution,  and  my  approximations,  to  the  Mechanical 
Equations.  Here  I  am  only  going  to  use  the  second  layer.  We  are  interested  in  the 
temperature  profile  between  the  2nd  and  3rd  layers  and  find  that  the  3rd  layer 
temperature  is  constant  with  the  bottom  of  the  2nd  layer  temperature.  I  still  need  the 
parameters  for  the  first  layer,  but  only  the  equations  for  the  2nd  layer.  The  reed 
parameters  are  the  same  as  above,  except  that  the  alphas,  gammas,  and  "little  d's"  are 
now  functions  of  L  and  io.  That  block  is  not  re-printed.  The  variable  functions  are 
now  re-written  as  functions  of  length  and  chopper  frequency. 
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■  The  Root  Solving  Block  -Approximations 

Recall,  above,  that  the  approximations  resulted  in  little  error  in  the  roots  and 
derivative  terms.  However,  one  trait  of  these  functions  is  that  a  slight  difference  in 
the  fit  can  change  the  answer  drastically.  Mathematica  calls  the  first  term  "  1 "  when 
it's  actually  "0"  so  to  calculate  the  proper  terms,  I  use  of  "k+1" 


lams[k_]  :=  |  (2  *  IntegerPart  [k]  +  1)  *  ; 

root  0  [k_]  :  =  -  2 . 02072  +  8 . 19207  *  IntegerPart  [k+1]  - 
8 .20216  *  (IntegerPart  [k  +  1] ) 2 ; 

.0001*  (root0[k] ) 


root[L_,  o) _ ,  k  ]  :  = 


w  +  L2 


sig[L_,  u_,  k_]  :=  sig[L,  co,  k]  =  -root[L,  u,  k]; 
derv[w_/  k_]  :  = 

(0.0009902715  -3.328841*  10‘9  *  IntegerPart  [  k]  - 

7 . 094576316  *  10~10  *  (IntegerPart [k] ) 2)  * 
(*the  derv  is  essentially  constant  with  L. *) 


wl[L_,  o)_,  k_]  :  = 

wl[L,  u,  k]  =  Ha[V  (root  [L,  co,  k]  +  ddl[L,  a>]  *  lams[k] )  ]  ; 
w2  [L_,  u_,  k_]  s  = 

w2[L,  Co,  k]  =  Re[V  (root  [L,  o),  k]  +  dd2  [L,  w]  *  lams[k] )  ]  ; 
w3  [L_,  o>_,  k_]  :  = 

w3[L#  u,  k]  =Re[V  (root [L,  <o,  k]  +dd3[L,  co]  *lams[k])]; 


The  Constant  Functions  "ck" 

ck2[L_,  (i) _ ,  k_]  :=  ck2  [L,  w,  k]  = 

2* 


ddl[L,  <u]  Sin[wl[L,  co,  k]  *gaml[L,  w]  ] 


dd2  [It,  <u] 


wl[L,  a,  k]  *derv[w,  k] 


■  The  Time  Behavior  "I" 

Incoming  Energy 
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kk=  1; 

(♦this  k  represents  an  absorption  of  the  top  layer*) 

Io=  1  *  10~6  (*  this  represents  the  initial  intensity 

of  the  incoming  beam,  I  figure  it1  s  reduced  quite 
a  bit  because  of  the  monochrometer  and  stuff  *)  ; 
to  =  . 5  (*  fraction  of  period  that  signal  is  on, 
due  to  chopper  setup  *) ; 

P[t_]  :=  P[t]  =  IntegerPart [t] 

Xi[L_,  w_,  t_,  k_]  :=Ii[L,  o>,  t,  k]  = 

kk*  Io  .  ri  . ,  _  I  (Esi«tL' _  l) 

_ *  ">*]  *t  ^  I  * 

sig[L,  u,  k]  \  w,k]  _  i) 

(EsifftL'  _  1)  _  E«iff[L,  j  + 

kk*  Xo 

-  - - -  *  ' 

sig[L,  a),  k] 

If[t-P[t]  <=  tO,  1,  E_Sig[L'  ]  ; 

■  The  X  Behavior,  <p,  or  "Thi" 

This  is  the  only  term  that  is  effected  by  integration  of  the  curvature  to  get  deflection. 

&(£) =  —-j  Sine  (Ik  £)  + 

Ak  Ak 


lam[k_]  :=  lam[k]  =  |  (2*  IntegerPart  [k]  +  1)  *  ; 

thiNew[f_,  k_]  :  = 

thiNew[f,  k]  =  - - - -  *Sin[lam[k]  *£]  + 

(lam[k] ) 3 
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■  The  Y  Behavior,  iff,  or  "Psi" 


Aa[L_,  w_,  k_]  :  = 

Aa[L,  0) ,  k]  =Cosh[w3[L,  u,  k]  *gam3[L,  w]  ]  ; 


Bb[L_,  0)_,  k_]  :  =  Bb[L,  w,  k]  = 
w3[L,  w ,  k] 


dd3  [L,  ft)] 


dd2[L,  ft)] 

*  Sinh[w3  [L,  ft>,  k]  *  gam3  [L,  ft)]  ]  ; 


w2[L,  ft),  k] 

Psi2 [rj_,  L_,  ft)_,  k_]  :=  Psi2[j7,  L,  w,  k]  =  Aa[L,  w,  k]  * 


Cosh[w2  [L,  to,  k]  * 
Bb[L,  0)g  k]  * 
Sinh[w2[L,  ft),  k]  * 


[  co] 


[  Vdd2[L/  co] 


+  gam2  [L,  oj] 


] 


+  gam2  [L,  ft)]  1  ]  ; 


■  Deflection 

Mathematica  is  very  strange  about  using  units.  It  does  not  like  to  cancel  units  even 
when  they  obviously  cancel,  so  occasionally  I  have  to  force  the  program  to  cancel, 
hence  use  of  such  terms  as  "  kelvin,  and  Meter2"  For  example,  Mathematica  will  not 

cance  Meti*  V  Meter2  .  The  deflection  as  a  function  of  frequency  is  adapted  from  the 
article  "Photoghermal  measurements  at  picowatt  resolution  using  uncooled 
micro-optomechanical  sensors"  Varesi,  Lai  et  al.  Notice  that  the  maximum  deflection 
does  NOT  occur  at  resonance  here.  I  expected  it  to  peak  there. 
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kmax  =  50; 

kb  =  Bolt  zmannConst  ant 


(*  1. 38065800000000017* 10~23  in  J/K  *) ; 

T  =  298  Kelvin  (*  roam  temperature*)  ; 

B=  8;  (*  Measurement  Bandwidth  *) 

Q  =  20;  (*  FWHM  of  peak  at  resonance  *) 
w 

kc [L_]  :  =  - *  (E3  *  t33  +  E2  *  t23) 

4*  (L  Meter) 3 

m  [L_]  :  = 

(♦effective*)  SI  [  .24* w*  (L  Meter)  *  (ro3  *  t3  +  ro2  *  t2)  ] 


wo[L_] 


kc[L] 

m[L] 


]* 


Second2 


(♦output  in  Hertz  *) 

zz[L_,  w_,  t_]  :=  zz [It,  w,  r]  =  6*  (alpha3  -  alpha2)  * 


(t3  + 12) 


kmax 


-I  *  /  ck2[L,  a),  k]  *I±[L,  (0,  z,  kl  * 
t2  *  K  )  Z_J 


k=0 


Psi2  [alphlfL]  # 
Meter  Kelvin  *  - 


L#  o)  , 


k]  *  thiNew[l#  k]  Meter 

=  *(WQ[L])2^ 

a>2  \ 

(WO[L])a  I 


(  (WO[L]  )  2  -  (CO)2)2* 


(wo[L])2*  (a)2  V 

0  J; 


■  Thermal  Noise 

RMS  Noise  (Thermally  Induced  Lever  Noise),  based  on  Sarid  Book. 


Wp[L_,  w_]  :=  Convert  [ 


4  *kb*T*B 
kc[L]  *mo[L]  *  Q 


(wo[L] ) 2  /  (y  (((wo[L]): 
"sj  Meter2  ]  * 


(w)2)  +  (wo[L] ) 2  *  (w)2)). 


Meter 
V  Meter2 


The  second  term  is  unitless  if  both  frequenceies  are  written  in  hertz,  so  the  unit  is 
removed  from  the  terms  that  already  have  it  to  appease  Mathematica.  Also,  the  terms 
below  are  canceled  in  terms  of  Meters,  assuring  that  SNR  is  the  correct  breed  of 
unitless. 
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SigNoise  [L_,  f_]  :=  SigNOise[L,  f]  = 


zz[L,  f,  20.5] 


Meter 


l  r  Q2  *kc[L]  *wo[L]  ,  , 

J Convert  f - -*  Meter2,  ll  ; 

V  1  4  *  kb*  T  J 


The  second  term  is  unitless  if  both  frequenceies  are  written  in  hertz,  so  the  unit  is 
removed  from  the  terms  that  already  have  it  to  make  MMA  happy. 


■  Sample  Output 

Plot  of  Deflection  vs  Chopper  frequency,  blow  up  to  show  that  the  peak  does  not 
occur  at  the  resonant  frequency  (shown  by  the  line). 

zz[  .0225,  f,  20.5] 

Plot  - , 

L  Meter 

{f,  1,  130},  GridLines  ->  { (wo[ .0225]  } ,  None}, 

PlotRange ->  {{95,  105},  Automatic}]; 


This  peaks  at  100,  with  a  max  at  3.1,  wo=  102.08 

.015  m  reed:  peaks  at  225.5,  at  again  about  3.25,  wo  =  229.7 

.01  m  reed:  peaks  at  508  at  3.25,  wo  =  516.8 

Difference  grows  inversly  proportional  to  length-  ~2, 4,  9  difference,  corresoponding 
to  a  2%,  1.8%  and  1.7%  drop. 

Now  to  generate  a  3D  plot  of  sensitivity: 
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Phase  Analysis 


■  Tl/2  points 

Here,  I  look  at  what  times  the  amplitude  is  one  half  of  maximum.  On  the  rise,  this  time  is 
called  te  and  on  the  fall  this  time  is  called  ti.  Phase  shift  will  be  recognized  by  a  change  in  the 
difference  between  these  two  times.  Using  the  time  varying  part  of  the  thermal  solution,  I,  the 
following  equation  is  obtained  for  that  difference  in  time,  as  a  function  of  the  root: 

tL-te  =  ^ln  (2+eT) 


delt[sig_]  :=  — —  *Log[2  +  Et2’1  ; 

sig  1  1 

Plot[N[delt[sig]  ] ,  {sig,  .001,  25}]; 


The  function  quickly  goes  to  the  constant  value  of  1/2.  Thus,  phase  information  can  be 
adequately  described  by  using  only  the  first  root  term.  There  is  some  loss  of  accuracy,  but  the 
general  trends  will  still  be  observed. 


■  Reed  Parameters  (same  as  in  previous  file,  not  reproduced) 


■  The  Root  Solving  Block  -  We  only  need  the  First  Root,  sig,  for  now. 

(same  as  previous  file,  only  root  value  reproduced  ) 

k  Ak  -d2Ak2  zk  Jj-(zk)  iUi(zk) 

0  -§-  -0.05608007  -0.05063511  +  0.1  0.002800636  +  0.1  0.  +  0.2183303 
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■  I  Analysis,  After  Many  Periods 

After  a  large  amount  of  periods  has  passed,  the  incoming  energy  deposition  function  can  be 
written  as: 

kk=l;  (*  absorption  *) 

Io=l;  (*  initial  intensity  *) 

to =  .5  (*  fraction  of  period  that  signal  is  on  *) ; 

P[t_]  :=  p[t]  =  IntegerPart [t] 

(★number  of  periods  that  has  passed  thus  far  *) 


(*  short  term  behavior  *) 
kk*  lo 

Ia[t_]  s= 
kk*  lo 


sig 


(/Eoig*to  _  \ 

— - —  *  (Esiff*p[t]  -  1)  -  Esia*p[t]  + 

(E«i9-1)  ) 


sig 


*If[t-p[t]  <=tO,  1,  E-8is,*((t-p^)-to)]  ; 


(*  long  term  behavior  *) 

kk*  lo  (  (Eslg*to  ■ 
Ib[t_]  :=  - - - *  1 


1) 


sig 


(Esig  -  1) 


*E‘“ig  >  + 


kk*  lo 
sig 


*If[t-p[t]  <=  to,  1,  E-^ttt-pttn-to)^ 


(★long  term  behavior,  not  modulated  by  the  constant  terms  *) 
lc[t_]  := 

((Esig*to  -  1)  \ 

(ESlg  _  -  Ij  * E-Si®  (t-pttn  +if[t-p[t]  <=  to,  l,  E-alsr*«t'p[tb-t°)  j . 

Plot[lc[t]  ,  {t,  0,  1}]; 

Plot[Ib[t]  ,  {t,  0,  1}]; 

Plot  [la  [t],  {t,  0,  1}]; 
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This  is  Strange...  the  expected  shape  is  not  evident.  Here  it  is  appears  linear! 

lb[ .5]  -  lb[0] 

0.249987 

The  amplitude  difference  in  one  period  is  seen  as  about  .25.  This  is  consistent  with  the 
following  analysis: 


■  I  Min  and  Max 

For  Large  Periods,  the  minimum  and  mazimum  values  of  the  time  dependant  portion  are  fairly 
constant.  The  most  important  quantity,  however  is  the  difference  between  the  two.  The  max, 
min,  and  the  difference,  as  a  function  of  the  root  itself,  can  be  obtained  from  the  long  term 
equation  as: 


j  -  _L 


i+«;t 

e"T 


1+e'T 


Amplitude  =  /max  - 1^  =  d-tanh(f-) 
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1 


1 


lmax[sig_]  :=  — ; - —  ; 

Slsr  i  +  b-t2- 


lmin[sig_] 
Amp[sig_]  :  = 


•la 


sig  i  +  e-¥- 

1  T  sig  .. 

- Tanh  — -1 

sig  1  4  1 


Plot [{Amp[ sig]  ,  Emaxfsig]  -Indn[sig]># 

{sig,  -2,  0},  AxesLabel  ->  {"root  value",  "Amplitude"}]; 
Plot [ {Anp [ sig] ,  Imax[sig]  -lmin[sig]},  {sig,  -1000,  0}, 
AxesLabel ->  {"root  value",  "Amplitude" } ] ; 


root  value 


The  Amplitude  starts  off  at  about  .25  for  the  smaller  roots,  and  quickly  approaches  zero  as  the 
roots  grow  in  magintude.  The  roots,  as  determined  previously,  start  off  at  a  small  negative 
number  and  then  progress  quickly  to  negative  numbers  far  beyond  the  range  of  the  plot.  Thus, 
using  the  first  root  is  a  valid  approximation  for  future  behaviour  analysis. 
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■  Resonance  and  Optimum  Thickness  Exploration 


■  Important  Results,  presented  upfront: 

wo[.0225  m]=  102.078  Hertz 
wo[.015  m]=  229.676  Hertz 
wo[.01  m]=  516.771  Hertz 

wo[.009  m]=  637.989  Hertz 


wo  in  Hertz 


Note  that  the  wider  the  reed  gets,  compared  to  the  same  length,  the  less  sensitive  the  reed  is 
overall. 


Sensitivity  — 


A1  Thickness 
Si  Thcickness 
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A  Close  Up  reveals  the  optimum  to  be  at  .325,  this  optimum  is  not  affected  by  length  or  width 
changes. 

Sensitivity  Change  with  length  and  modulation  frequency 


Sensitivity 


m 

¥ 


As  length  increases,  resonant  frequency  increases  and  sensitivity  decreases.  If  the  same  amount 
of  energy  is  absorped  in  the  reed  for  all  three  lengths  (.02  m,  .015  m,  and  .01  m)  then  the 
deflection  behavior  will  be  similar  to  this  sensitivity  relationship 
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This  3-D  graphs  shows  the  effect  of  chopper  frequency  and  length  of  the  reed  on  the 
sensitivity .  Notice,  from  the  shape  of  the  graph,  that  current  reeds  in  testing  have  resonant 
frequencies  so  low  that  they  are  already  operating  in  their  resonant  frequency,  their  most 
sensitive  area! 


■  Reed  Parameters 

<<  Miscellaneous' units' ; 

«  Miscellaneous ' PhysicalConstants ' 

L  =  .02  Meter; 
w  =  .  002  Meter; 


(*  Layer  one,  M&A1204:  Hartnett’ s  values  *) 
tl  =  1  *  10~6  Meter; 


rol  =  3.58 


Gram 


El  =  2.77  *  1011 


Centimeter3 
Newton 


Meter2 


(*elastic  modulus*) 
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Watt 

laml  =14,6 - (*  thermal  conductivity*)  ; 

Meter  *  Kelvin 

.  1 

alphal  =6.5*  10'  - : — ;  (* thermal  expansion  for  30-00C*) 

Kelvin 

Calorie 

cpl  =  .200  - ;  (*at  20  C  Roy1  s  value*) 

Gram*  Celsius 


.  _  r  laml  Micro  Meter2  - 

Dif f  1  =  Convert  - ,  -  ; 

L  cpl  *  rol  Second  J 


(*  Layer  two  Silicon  *) 


t2  =  8  *  10'6  Meter; 

Kilogram 

ro2  =  2328  - , 

Meter3 
--  Newton 


E2  =  1  *  1011 


Meter2 


lam2  =  150 


alpha2  =  2.6*  10 


Meter  Kelvin 
_  1 


cp2  =  700 


Kelvin 

Joule 

Kilogram  Kelvin 


Dif £2  =  N [Convert [ 


lam2  Micro  Meter2 
cp2  *  ro2  9  Second 


(*  Layer  three  Aluminum  *) 
t3  =  .  5  *  10"6  Meter; 

Kilogram 


ro3 =  2702 


E3  =  ^lO1 


lam3  =  237 


alpha3  =  23.6*10 


Meter3 
--  Newton 

1  Meter2  7 
watt 

Meter  Kelvin 

_  -  1 


(*  thickness*) 
(*  density  *) 


(*  Elastic  Modulus  *) 


(* thermal  conductivity*) 


-;  (* thermal  expansion*) 


Kelvin 

*  Joule 

cp3  =  908  — — - : — ;  (*heat  capacity*) 

Kilogram  Kelvin 


Diff3 = N [convert [ 


lam3  Micro  Meter2 
cp3  *  ro3  9  Second 


Printout  of  the  Diffusion  Coefficients: 
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Diffl 

Di££2 

Dif£3 

4 . 87032  Meter2  Micro 
Second 

92 . 0471  Meter2  Micro 
Second 

96 . 6  Meter2  Micro 
Second 


■  Mechanical  Resonance  Prediction 

w 

n_]  s=  - *t23  (E3  *  n3  +  E2) 

4*  (L  Meter) 3 

m  [L_#  n_]  :  =  (^effective*) Si [ .24 *  w*  (LMeter)  * t2  (ro3*n+ro2)] 

wo[L_,  n_]  :=  Convert f  / Sl[ —  L'  1  *  Second2  #  l]  Hertz; 

L  V  Lm[L,  n]  J  J 

wo[ .0225 ,  .0625] 

102 . 078  Hertz 


Plot[ 


wo  [ .0225,  n] 
Hertz 


{n#  0,  1}, 


Al  Thickness 

_ _  n 

Si  Thickness 
AxesOrigin  ->  {0,  87}]; 


"wo  in  Hertz"}, 


wo  in  Hertz 
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-  wo[L,  n] 

Plot3D[ - ,  {L,  .01,  . 03}  ,  {n,  0,  1},  PlotPoints  ->  50l  ; 

Hertz  J 


Be  sure  to  note  that  n  is  the  ratio,  (n  =  )  The  change  of  Resonant  frequency 

WRT  length  has  the  same  shape  for  each  thickness  of  Aluminum  layer.  As  is  better  seen  in  the 
2D  graph,  the  resonant  frequency  hits  a  minimum  at  about  n=  .6.  This  means  for  minimum 
resonant  frequency,  we'd  want  to  operate  with  an  aluminum  layer  that  was  4.8  microns  thick. 
We  may  want  to  operate  at  minimum  resonant  frequency,  because  the  tests  show  that  smaller 
reeds  are  more  stable-  also,  the  equipment  may  be  more  stable  if  we  can  bring  the  resonant 
frequency  down. 


■  Optimum  Thickness  Ratio  Prediction 

■  thickness  ratio  (n)  effect  on  sensitivity: 

Here  the  effect  on  varying  the  thickness  of  the  Aluminum  layer  with  respect  to  that  of  the 
Silicon  layer  is  explored.  The  silicon  layer  is  held  constant  because  it  is  the  pre-purchased 
material.  The  ratio  n  =  £or  4^:ic,kness . 

t2  Si  thickness 
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(*  Young1  s  Mbdulus  *) 
tl 


n  = 


t2 


El 

thi  =  - ; 

E2 


K=4  +  6*n+4*n2  +  thi  *  n3  + 


thi  *n 


7.706 


(*  Several  Other  Ratios  of  that  are  handy:  *) 


gamma  = 

beta  = 


lam3 


lam2 

alpha3 


theta  = 


alpha2  ' 

(beta  -  1)  *  (n+  1) 


K*  (gamma* n+  1) 


psi  = 


ro3 


ro2 

S[n_,  L_,  w_]  s  = 


alpha2  L3 
2  * -  - * 


(beta  -  1)  *  (n  +  1) 


lam2  t22  *w  (4  +  6*n+4*n2  +  thi*n3  +  )  *  (gamma*  n+  1) 

Watt 


„  r  /  Watt  \ 

Plotl  S[n,  L,  w]  * -  ,  {n,  .01,  1},  PlotPoints  ->  100, 

v  Meter  / 


Meter 

_  „  f  A1  Thickness 
AxesLabel  ->  { "n=  — 


",  "Sensitivity  ■—■n}] 


Si  tfhcickness 


W 


Sensitivity  — 


A1  Thickness 
Si  Thcickness 


A  close  up  reveals  the  optimum  to  be  at  0..325: 
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■  Sensitivity  Optimum  Other  Ratios?  Conductivity,  Expansion? 

S[n_,  ga_]  :  = 

alpha2  L3  (beta  -  1)  *  (n  +  1) 

2  *  -  - *  — - — - • 

lam2  t22*w  (4  +  6*n+ 4*n2  +  thi*n3  + *  (ga*n+ 1) 

Watt 

a  =  Plot [S[ .0325,  ga]  * - ,  {ga,  .01,  20}  (*,  PlotPoints->100*) 

Meter 

(*,AxesLabel->{"ga=  4r^^-n, "Sensitivity  £"}*)], 

,  „  r  Watt 

b  =  Plot  I S[ .325,  ga]  * - ,  {ga,  .01,  20}  (*,  PlotPoints->100*) 

L  Meter 

(*,AxesLabel->{ "ga=  "."Sensitivity  £"}*)]; 

,  r  Watt 

c =  Plot [S[ .5,  ga]  * - ,  {ga,  .01,  20}  (*,  PlotPoints->100*) 

Meter 

(*,AxesLabel->{nga=  " , "Sensitivity  £"}*)]; 

.  Watt 

d  =  Plot [S [1,  ga]  * - ,  {ga,  .01,  20}  (*,  PlotPoints->100*) 

L  Meter 

( * , AxesLabel- > { " ga=  jf  g°^°;  " ,  "Sensitivity  £»}*)]; 

Show  [a,  b,  c,  d]  ; 

Below  is  a  plot  of  sensitivity  variance  with  conductivity  ratios,  the  conductivity  of  the  bottom 
material  divided  by  the  silicon  conductivity. 
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Plot  [s[n,  1]  * - ,  {n,  .01,  1}  1  ; 

Meter  J 


Plot  [S[n,  1]  * - ,  {n,  .01,  1}, 

Meter 

PlotRange  ->  { { .2,  .45},  Automatic}]; 


For  a  conductivity  ratio  of  1,  the  optimum  thickness  ratio  changes  to  n  =  .367.  Also,  went  up 
from  1.495  at  max  to  about  1.718  This  makes  sense.  If  you  make  the  bottom  layer  less 
conductive,  it  will  require  more  of  itself  to  keep  in  balance. 


■  Sensitivity  as  a  function  of  frequency  (assuming  w  =  .1  L)  for  three  different  lengths 

Also  Assuming  Q=20,  a  value  that  seems  reasonable  judging  from  the  literature.  This  must 
actually  be  measured  in  the  lab  once  the  resonance  peak  has  been  graphed  out  and  is  dependent 
on  resonant  frequendcy,  thus  also  length  of  the  reed  though  I've  assumed  it  constant. 

Reed  Parameters  block,  identical  to  that  printed  earlier  in  this  file,  has  been  hidden  for  printing. 
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■  2  cm  Reed  (the  equations  are  reproduced  only  for  this  case,  but  similar  ones  were  used 
to  produce  the  other  peaks  on  the  graph) 

Units  are  thrown  around  and  balanced  because  Mathematica  does  not  like  to  operate  on  units  as 
well  as  it  could. 


L  =  .02  Meter; 
w  =  .1  *  L; 

k  =  - - - *  (E3  *  t33  +  E2  *  t23) 

4*  (L)3 

m=  (*effective*)SI[.24*w*  (L)  *  (ro3*t3  +  ro2*t2)  ] 
k 


m 


S  [om_]  :  =  2  *  (alphal  -  alpha2)  * 
L3 


I  (tl  +  t2)  \ 

i  t22  *  K  i 


wo' 


(laml*  tl  +  lam2  *  t2)  *w  _^(wo2  _  ^  2  + 
Sen [4  Hertz] 


wo^om2 

Q2 


Hertz4 


Hertz2 


0.00320062  Newton 
Meter 

1 . 9176  x  10"7  Kilogram 

16690 . 8  Newton 
Kilogram  Meter 

129 . 193  Hertz 

0 . 469214  Meter 
Watt 
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Show[a,  b,  c,  AxesLabel  ->  {"Mod.  Freq.  in  Hz",  "Sensitivity 


Sensitivity  — 
W 


■  3D  Plot,  length  along  one  axis,  and  frequency  along  another,  assuming  w  =  .1  L. 

The  first  half  of  the  equation  is  normal,  but  the  unit  play  in  the  second  half  is  because  of 
Mathematical  intrinsic  funnyness.  Originally,  this  program  was  chosen  for  ease  of  use,  but 
this  is  only  the  case  for  scaled  work. 


w[L_]  :=  .  1*L; 

w[L] 

kk[L_]  :=  - —  *  (E3  *  t33  +  E2  *  t23) 

4  *  (L  Meter ) 3 

rnn[L_]  :=  (*effective*)  SI[  .24  *  w[L]  *  (LMeter)  *  (ro3  *  t3  +  ro2  *  t2)  ] 

wwo[L_]  :=  Convert  L/  sif-^— L]  *  Second2  ,  ll  Hertz; 

LV  L  irm[L  ]  1  1 


Sen[am_,  L_]  s  = 


2  *  (alphal  -  alpha2)  * 


(tl  +  t2)  \  L3 

-  * - - —  * 

t22  *  K  /  (laml*  tl  +  lam2  *  t2)  *w[L] 


Sen [ 4  Hertz,  .  02  Meter] 

0 . 469214  Meter 
Watt 
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Plot 3D [ 

,  Watt 

Sen[om Hertz,  L Meter]  * - ,  {can,  1,  600),  {L,  .005,  .025}, 

Meter 


Axes  Label  ->  ("mod.  freq.  Hz",  "Length,  m",  "Sensitivity  - 

nW  J 

PlotPoints  ->  50] ; 


C.43 


■  More  Resonance  Plots 


<<  Graphics 'Legend'; 

k[L_,  n_]  :=  - - - *t23  (E3*n3  +  E2) 

4  *  (L  Meter ) 3 

m  [L_»  n_]  :=  (*effective*)  SI[  .24  *w*  (LMeter)  *t2  (ro3*n+ro2)] 


wo[L_,  n_]  ;=  Convert  L/  Slf-^L'  —  -  1  *  Second2  ,  ll  Hertz; 

LV  Lm[L,  n]J  J 

pi.t[{^1L_-  •°625) , 


wo[L, 


Hertz 
25]  vto[L, 


.6]  wo[L,  1] 


},  {L,  .009,  .0225}, 


Hertz  Hertz  Hertz 

Frame  ->  True,  FrameLabel  ->  {nL  (m)  ",  "wo  in  Hertz"}, 

AxesOrigin  ->  {  .  01,  200}  ,  PlotStyle  ->  {Dashing[  {  .1,  .  1}  ]  , 

Dashing  [  {  .05,  .05}],  Dashing[{  .01,  .05}],  RGBColor[0,  0,  0],}, 
PlotLegend  ->  {"n  =  .0625",  "n  =  .25",  "n  =  .6",  "n  =  1”}, 
LegendPosition  ->  {1,  - .  5}  , 

LegendShadow ->  None,  LegendSize ->  {.5,  1}, 

PlotRange  ->  {{.01,  .015},  {200,  530}},  Plot Label -> 

"Resonant  Frequency  vs  Length  (in  efficient  length  range) ■  ] ; 


Resonant  Frequency  vs  Length  (in  efficient  length  range) 


L  (m) 


■  Estimate  for  Deflection  due  to  constant  1.269*  10“6  amplitude  temperature  vibrations. 

Using  lab  data  to  approximate  Io,  temperature  fluctuations  between  1 .269  *  10~6  Kelvin, 
and  1.295*  1(T6  Kelvin  were  calculated  forthe  2.25  cmfour- 
layer  reed .  Below  is  the  corresponding  deflection : 


C.44 


(t2  +  t3) 


*  L2  *  delT; 


z  [delT_]  :  =  3  *  (alpha3  -  alpha2)  * 

t22  *  K 

delt  =  1.295  *  1<T6  Kelvin  -  1.269  *  1<T6  Kelvin 
Convert  [ z  [delt] ,  Nano  Meter] 

2 . 6  x  1CT8  Kelvin 

0 . 0112923  Meter  Nano 

This  is  small  compared  to  the  ~  1-10  nm  measurements  currently  observed  in  the  lab. 
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